CN102520726B - Estimation method of atmosphere angle of attack and angle of sideslip in high-angle-of-attack flight status - Google Patents

Abstract

The invention discloses an estimation method of an atmosphere angle of attack and an angle of sideslip in a high-angle-of-attack flight status. Under the premise that the structure of an aircraft is not changed and additional measuring devices and hardware devices are not added, the parameters output by the existing onboard strapdown inertial navigation system, a GPS (Global Position System) system, a flight control system and a fuel measuring system are fully used for resolving the angular acceleration, the rotational inertia, the product of inertia, the coefficient of a moment equation set and the total mass of the aircraft; based on a flight dynamics model, by establishing an extended Kalman filter, and integrating the processes of resolving the flight dynamics model and the status estimation together, accurate estimation of the current angle of attack, the angle of sideslip and the true airspeed is realized, and then the current air speed and the change rate of the air speed are resolved and input into the extended Kalman filter as the feedback information so as to finish the real-time accurate estimation of the angle of attack and the angle of sideslip of next time; the real-time accurate estimation of the angle of attack and the angle of sideslip is realized by a recursive resolving mode, measurement range and measurement precision of an atmosphere data system are improved, and adaptability to complex flight environment is also improved.

Description

Atmospheric attack angle and sideslip angle estimation method in large attack angle flight state

Technical Field

The invention relates to an atmospheric attack angle and sideslip angle estimation method in a large attack angle flight state, and belongs to the technical fields of atmospheric data, aircraft control, inertial navigation, GPS (global positioning system) and the like.

Background

The attack angle and the sideslip angle are two important parameters of the flight state and the flight mechanics of the aircraft, and the precision and the reliability of the attack angle and the sideslip angle are directly related to the normal work and the performance exertion of systems of aircraft control, navigation, an engine, fire control and the like. The main tactical performance requirements of modern high performance fighters include ultra-maneuvering capability, supersonic cruise, stealth, high information advantage, etc. The high maneuvering flight means that the aircraft has the flight capability of a large attack angle, the airframe surface airflows are separated due to the flight of the large attack angle, the aircraft is in the nonlinear vortex airflows, the current attack angle and sideslip angle are mainly measured through a vane sensor, a differential pressure sensor, a zero differential pressure sensor and the like which are installed on the aircraft, the attack angle and sideslip angle sensors are sensitive to airflow disturbance, the performances of the attack angle and sideslip angle sensors are seriously reduced, the accurate measurement of the attack angle and the sideslip angle in the flight range of the large attack angle is difficult to meet, and the performance exertion of the high-performance aircraft is seriously restricted. In addition, the outstanding sensor measuring device will reduce the stealth of the aircraft, so the problem of measuring the attack angle and the sideslip angle of the high-performance aircraft is very important to solve.

At present, research on measuring methods of a large-attack-angle lower attack angle and a large-sideslip angle is carried out at home and abroad. In 1991, a 'Development of a Pneumatic High-Angle-of-attach flash air description Sensing (HI-FADS) System' published by NASA report of America, a built-in atmosphere data System for measuring the pressure distribution of the surface of an aircraft by means of pressure sensor arrays embedded at different positions at the front end of the aircraft and estimating atmospheric parameters from the pressure distribution is proposed, and the method can realize the measurement of an Attack Angle and a sideslip Angle under a large Attack Angle flight environment; 2009 us patent US7,564,539B2: the Optical Air Data Systems Methods provide an Optical atmospheric Data system which adopts laser, atmospheric molecules and aerosol to respectively generate Rayleigh scattering and Mie scattering effects and calculates atmospheric parameters according to the strength of a backscattering signal and Doppler frequency shift, and the method is not influenced by a flight state when measuring the atmospheric parameters and meets the measurement of an attack angle and a sideslip angle under a flight state with a large attack angle; the two methods for researching the novel embedded atmospheric data system and the optical atmospheric data system are high in technical difficulty and high in cost; in 2006, an article "Flight Testing of the X-45AJ-UCAS Computational Alpha-Beta System" published by AIAA in the Conference "AIAA guiding, Navigation, and Control Conference and inhibition" provides a method for estimating an attack angle and a sideslip angle by using an inertial Navigation System, a Flight Control System and an atmospheric data System, so that the influence of the characteristics of a large attack angle on the measurement of the Flight attack angle and the sideslip angle is avoided, but the method does not consider the influence of the actual wind speed around an aircraft, and the estimated attack angle and sideslip angle have System errors.

Disclosure of Invention

The invention aims to provide an atmospheric attack angle and sideslip angle estimation method under a large attack angle flight state by fully utilizing output parameters of an airborne strapdown inertial navigation system, a GPS (global positioning system), a flight control system and a fuel oil measurement system and based on a flight dynamics model aiming at the defects of the existing large attack angle atmospheric parameter measurement technology.

In order to achieve the purpose, the invention adopts the following technical scheme:

an atmospheric attack angle and sideslip angle estimation method in a large attack angle flight state is characterized by comprising the following steps: parameters output using equipment onboard an aircraft, including:

(1) attitude angle, angular velocity and linear acceleration output by the strapdown inertial navigation system;

(2) ground speed output by the GPS system;

(3) aerodynamic force, external torque and engine thrust output by the flight control system;

(4) the mass of the residual fuel output by the fuel measuring system;

solving angular acceleration through a differential method according to the angular velocity;

solving the rotational inertia and the inertia product of the aircraft by the angular velocity, the angular acceleration and the external moment through a least square estimation method, and further obtaining the coefficient of a moment equation set;

adding the mass of the residual fuel oil with the known mass of an aircraft body, the mass of airborne equipment, the mass of passengers and the mass of weapons to obtain the total mass of the aircraft;

constructing an extended Kalman filter according to the flight dynamics model, and realizing accurate estimation of an attack angle, a sideslip angle and a vacuum speed at the current moment in a large attack angle flight state by using an extended Kalman filtering algorithm;

calculating the wind speed at the current moment according to the attitude angle and the ground speed at the current moment and the estimated attack angle, the sideslip angle and the vacuum speed at the current moment, and calculating the wind speed change rate at the current moment by combining the wind speed at the previous moment;

and feeding the obtained wind speed and the wind speed change rate at the current moment back to the extended Kalman filter mode algorithm module for finishing the real-time accurate estimation of the attack angle and the sideslip angle at the next moment, and realizing the real-time accurate estimation of the attack angle and the sideslip angle in a recursion solving mode.

The method comprises the following steps:

(1) reading three angular velocities, three attitude angles and three linear acceleration information in a strapdown inertial navigation system by using a period delta T, wherein the three angular velocity information are roll angular velocity p, pitch angular velocity q and yaw angular velocity r in the directions of an x axis, a y axis and a z axis under a machine body coordinate system respectively; the three attitude angle information are a pitch angle theta, a roll angle phi and a yaw angle psi respectively; the three pieces of linear acceleration information are respectively the linear acceleration a in the x-axis direction under the body coordinate system_xLinear acceleration a in the y-axis direction_yLinear acceleration a in the z-axis direction_zThe directions of an x axis, a y axis and a z axis of the body coordinate system are forward, rightward and downward respectively;

(2) reading the ground speed under the geographic coordinate system output by the GPS system by a period delta T, wherein the directions of the geographic coordinate system are defined as north, east and ground, and the ground speeds in the three directions are u respectively_g、v_g、w_g；

(3) Reading the aerodynamic force, the external torque and the engine thrust of the aircraft output by the flight control system at a period delta T, wherein the aerodynamic force of the aircraft is x under an airflow coordinate system_aAxial negative direction, y_aAxial direction and z_aThe components in the negative direction of the shaft are respectively resistance D, lateral force Y and lift L, and the components in the directions of the x-axis, the Y-axis and the z-axis of the thrust of the engine under a body coordinate system are respectively T_x、T_y、T_zThe components of the external torque in the directions of the x axis, the y axis and the z axis under the body coordinate system are rolling torque respectively

Pitching moment M and yawing moment N;

(4) reading the residual fuel mass m output by the fuel measurement system of the aircraft in a period delta T_f；

(5) According to the roll angular velocity p, the pitch angular velocity q and the yaw angular velocity r under the body coordinate system obtained in the step (1), wherein three angular velocity information at the time T are the roll angular velocity p (T), the pitch angular velocity q (T) and the yaw angular velocity r (T), three angular velocity information at the time T + delta T are the roll angular velocity p (T + delta T), the pitch angular velocity q (T + delta T), the yaw angular velocity r (T + delta T), three angular acceleration information in the directions of the x axis, the y axis and the z axis under the body coordinate system at the time T, and the roll angular acceleration

Acceleration of pitch angleYaw angular acceleration

Calculated by the differential method, namely:

$\left\{\begin{array}{c}\stackrel{\·}{p}\left(t\right)=\frac{p(t+\mathrm{\ΔT})-p\left(t\right)}{\mathrm{\ΔT}}\\ \stackrel{\·}{q}\left(t\right)=\frac{q(t+\mathrm{\ΔT})-q\left(t\right)}{\mathrm{\ΔT}}\\ \stackrel{\·}{r}\left(t\right)=\frac{r(t+\mathrm{\ΔT})-r\left(t\right)}{\mathrm{\ΔT}}\end{array}\right.;$

(6) known aircraft angular motion equation set $\left\{\begin{array}{c}\stackrel{\·}{p}{I}_x-\stackrel{\·}{r}{I}_{\mathrm{xz}}+\mathrm{qr}(I_z-I_y)-\mathrm{pq}{I}_{\mathrm{xz}}=\stackrel{\‾}{L}\\ \stackrel{\·}{q}{I}_y+\mathrm{pr}(I_x-I_z)+(p^2-r^2)I_{\mathrm{xz}}=M\\ \stackrel{\·}{r}{I}_z-\stackrel{\·}{p}{I}_{\mathrm{xz}}+\mathrm{pq}(I_y-I_x)+\mathrm{qr}{I}_{\mathrm{xz}}=N\end{array}\right.$ Is added to obtain the equation EJ ═ U, where $E=[\stackrel{\·}{p}+\mathrm{pr}-\mathrm{pq},\stackrel{\·}{q}-\mathrm{qr}+\mathrm{pq},\stackrel{\·}{r}+\mathrm{qr}-\mathrm{pr},p^2-r^2-\stackrel{\·}{p}-\stackrel{\·}{r}-\mathrm{pq}+\mathrm{qr}],$ J＝[I_x，I_y，I_z，I_xz]^T′、

T' denotes transpose, I_x、I_y、I_zThe inertia products of the aircraft to the x axis and the y axis are I_xyThe product of inertia about the x-axis and z-axis is I_xzThe product of inertia about the y-axis and the z-axis is I_yzThe aircraft is symmetrical about the plane of symmetry Oxz of the coordinate system of the aircraft body, the product of inertia I_xy＝I_yz＝0，I_xzIs non-zero; from t_aTime begins to t_bIn the time period when the moment is over, continuously acquiring n groups of roll angular velocities p, pitch angular velocities q and yaw angular velocities r by using the step (1), and continuously acquiring n groups of roll torques by using the step (3)

Pitching moment M and yawing moment N, and continuously obtaining N groups of roll angular accelerations in step (5)

Acceleration of pitch angleYaw angular acceleration

Obtaining n groups of roll angular velocity p, pitch angular velocity q, yaw angular velocity r and roll angular acceleration

Acceleration of pitch angle

Yaw angular acceleration

Rolling torque

Pitching moment M, yawing moment N and equation EJ ═U, using least square estimation method to calculate the moment of inertia I of the aircraft_x、I_y、I_zProduct of sum and inertia I_xzWherein n > 4;

(7) the moment of inertia I obtained according to the step (6)_x、I_y、I_zProduct of sum and inertia I_xzFor known aircraft moment equations $\left\{\begin{array}{c}\stackrel{\&CenterDot;}{p}=(c_{1}r+c_{2}p)q+c_3\stackrel{\&OverBar;}{L}+c_{4}N\\ \stackrel{\&CenterDot;}{q}={\mathrm{<figure-callout\; id="5"\; label="c"\; filenames="HDA0000121857970000032.png"\; state="\{\{state\}\}">c</figure-callout>}}_5\mathrm{pr}-c_6(p^2-r^2)+c_{7}M\\ \stackrel{\&CenterDot;}{r}=(c_{8}p-c_{2}r)q+c_4\stackrel{\&OverBar;}{L}+c_{9}N\end{array}\right.$ 9 coefficients c in₁、c₂、c₃、c₄、c₅、c₆、c₇、c₈、c₉Carry out a solution in which ${\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="HDA0000121857970000012.png,HDA0000121857970000041.png"\; state="\{\{state\}\}">c</figure-callout>}}_1=\frac{(I_y-I_z)I_z-I_{\mathrm{xz}}^2}{\mathrm{\&Sigma;}},$ $c_2=\frac{(I_x-I_y+I_z)I_{\mathrm{xz}}}{\mathrm{\&Sigma;}},$ $c_3=\frac{I_z}{\mathrm{\&Sigma;}},$ $c_4=\frac{I_{\mathrm{xz}}}{\mathrm{\&Sigma;}},$ ${\mathrm{<figure-callout\; id="5"\; label="c"\; filenames="HDA0000121857970000032.png"\; state="\{\{state\}\}">c</figure-callout>}}_5=\frac{I_z-I_x}{I_y},$ $c_6=\frac{I_{\mathrm{xz}}}{I_y},$ $c_7=\frac{1}{I_y},$ $c_8=\frac{I_x(I_x-I_y)+I_{\mathrm{xz}}^2}{\mathrm{\&Sigma;}},$ $c_9=\frac{I_x}{\mathrm{\&Sigma;}},$ Wherein $\mathrm{\&Sigma;}=I_x{I}_z-I_{\mathrm{xz}}^2;$

(8) The mass m of the residual fuel oil obtained in the step (4)_fAdding the mass of the aircraft body, the mass of airborne equipment, the mass of passengers and the mass of weapons to obtain the total mass m of the aircraft;

(9) selecting a pitch angle, a roll angle speed, a pitch angle speed, a yaw angle speed, an attack angle, a sideslip angle and a vacuum speed of the aircraft as state quantities according to the flight dynamics model, and further establishing an extended Kalman filter state equation; selecting a pitch angle, a roll angular velocity, a pitch angle velocity, a yaw angular velocity, a three-dimensional linear acceleration and a vacuum velocity as measurement quantities, and further establishing an extended Kalman filter measurement equation; according toThe current time t acquired in the step (1) and the step (2)_k+1Measuring information of time, the last time t obtained in step (3)_kAerodynamic force, external torque and engine thrust at the moment, 9 torque equation set coefficients obtained in step (7) and t obtained in step (8)_kTotal mass of the aircraft at time, and t_kThe three-dimensional wind speed and the three-dimensional wind speed change rate at the moment are obtained by using an extended Kalman filter equation_k+1The optimal estimated value of the state quantity at the moment is realized, so that t in the flight state with a large attack angle is realized_k+1Real-time accurate estimation of the angle of attack and sideslip angle at a time, where t_k＞t_b；

(10) Obtaining t according to step (1)_k+1Pitch angle, roll angle, yaw angle at the moment and t obtained in step (2)_k+1Combining the three-dimensional ground speed parameter of the moment with the t estimated in the step (9)_k+1The attack angle, the sideslip angle and the vacuum speed at the moment are calculated to obtain t_k+1Three dimensional wind speed at time, combined with t_kThe three-dimensional wind speed at the moment is calculated to obtain t_k+1The three-dimensional wind speed rate of change at the moment;

(11) t obtained in the step (10)_k+1Feeding back the three-dimensional wind speed and the three-dimensional wind speed change rate at the moment to the extended Kalman filtering module algorithm module for completing the next moment, namely t in the step (9)_k+2And estimating the moment attack angle and the sideslip angle.

The Kalman filtering algorithm in the step (9) comprises the following specific steps:

(a) establishment of extended Kalman filter state equation

According to the flight dynamics model, the pitch angle theta, the roll angle phi, the roll angle speed p, the pitch angle speed q, the yaw angle speed r, the attack angle alpha, the sideslip angle beta and the vacuum speed V of the aircraft are selected as state quantities X, namely X is [ theta, phi, p, q, r, alpha, beta, V]^T′And further establishing a state equation of the extended Kalman filter:

$\stackrel{\&CenterDot;}{X}\left(t\right)=f[X\left(t\right),u\left(t\right),t]+G\left(t\right)w\left(t\right)$

the specific form of the equation of state is:

wherein u (t) is control input quantity, and three-dimensional wind speeds along x-axis, y-axis and z-axis directions under a body coordinate system are respectively set as u_w、v_w、w_wThree-dimensional wind speed u_w、v_w、w_wRespectively isSystem noise $w\left(t\right)={{[\mathrm{\&delta;}}_{\stackrel{\&OverBar;}{L}},{\mathrm{\&delta;}}_M,{\mathrm{\&delta;}}_N,{\mathrm{\&delta;}}_D,{\mathrm{\&delta;}}_Y,{\mathrm{\&delta;}}_L,{\mathrm{\&delta;}}_{T_x},{\mathrm{\&delta;}}_{T_y},{\mathrm{\&delta;}}_{T_z},{\mathrm{\&delta;}}_{u_w},{\mathrm{\&delta;}}_{v_w},{\mathrm{\&delta;}}_{w_w},{\mathrm{\&delta;}}_{{\stackrel{\&CenterDot;}{u}}_w},{\mathrm{\&delta;}}_{{\stackrel{\&CenterDot;}{v}}_w},{\mathrm{\&delta;}}_{{\stackrel{\&CenterDot;}{w}}_w},{\mathrm{\&delta;}}_m]}^{T^{\&prime;}}$ Respectively by external torque of aircraft

M, N, aerodynamic force D, Y, L, engine thrust T_x、T_y、T_zWind speed u_w、v_w、w_wRate of change of wind speed

And total mass m, g is the acceleration of gravity, and the system noise figure array g (t) is defined as:

$G\left(t\right)=\left[\begin{array}{cccc}{0}_{2\&times;3}& 0_{2\&times;6}& 0_{2\&times;3}& 0_{2\&times;4}\\ A_{3\&times;3}& 0_{3\&times;6}& 0_{3\&times;3}& 0_{3\&times;4}\\ 0_{3\&times;3}& B_{3\&times;6}& C_{3\&times;3}& D_{3\&times;4}\end{array}\right]$

wherein, $A_{3\&times;3}=\left[\begin{array}{ccc}{c}_3& 0& c_4\\ 0& c_7& 0\\ c_4& 0& c_9\end{array}\right]$

$B_{3\&times;6}=\left[\begin{array}{cccccc}0& 0& \frac{-1}{\mathrm{mV}\cos \mathrm{\&beta;}}& \frac{-\sin \mathrm{\&alpha;}}{\mathrm{mV}\cos \mathrm{\&beta;}}& 0& \frac{\cos \mathrm{\&alpha;}}{\mathrm{mV}\cos \mathrm{\&beta;}}\\ 0& \frac{1}{\mathrm{mV}}& 0& \frac{-\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{\mathrm{mV}}& \frac{\cos \mathrm{\&beta;}}{\mathrm{mV}}& \frac{-\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{\mathrm{mV}}\\ \frac{-1}{\mathrm{<figure-callout\; id="0"\; label="m"\; filenames="HDA0000121857970000031.png,HDA0000121857970000032.png"\; state="\{\{state\}\}">m</figure-callout>}}& 0& 0& \frac{\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}}{m}& \frac{\sin \mathrm{\&beta;}}{m}& \frac{\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}}{m}\end{array}\right]$

$C_{3\&times;3}=\left[\begin{array}{ccc}\frac{q\cos \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}& \frac{-p\cos \mathrm{\&alpha;}-r\sin \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}& \frac{q\sin \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}\\ \frac{-r\cos \mathrm{\&beta;}-q\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}& \frac{-r\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}+p\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}& \frac{p\cos \mathrm{\&beta;}+q\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}\\ -r\sin \mathrm{\&beta;}+q\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}& r\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}-p\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}& p\sin \mathrm{\&beta;}-q\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}\end{array}\right]$

$D_{3\&times;4}=\left[\begin{array}{cccc}\frac{\sin \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}& 0& \frac{-\cos \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}& \frac{L+T_x\sin \mathrm{\&alpha;}-T_z\cos \mathrm{\&alpha;}}{m^{2}V\cos \mathrm{\&beta;}}\\ \frac{\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}& \frac{-\cos \mathrm{\&beta;}}{V}& \frac{\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}& \frac{-Y+T_x\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}-T_y\cos \mathrm{\&beta;}+T_z\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{m^{2}V}\\ -\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}& -\sin \mathrm{\&beta;}& -\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}& \frac{D-T_x\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}-T_y\sin \mathrm{\&beta;}-T_z\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}}{m^2}\end{array}\right]$

(b) establishment of extended Kalman filter measurement equation

Three-dimensional ground speed u obtained from GPS_g、v_g、w_gBy using

The calculated V is used as the vacuum speed measurement and is combined with the pitch angle theta, the roll angle phi, the roll angle speed p, the pitch angle speed q, the yaw angle speed r and the three-dimensional linear acceleration a which are obtained by the strapdown inertial navigation system_x、a_y、a_zThe quantity of the total composition Z, i.e. Z ═ θ, φ, p, q, r, a_x，a_y，a_z，V]^T′And further establishing an extended Kalman filter measurement equation:

Z(t)＝h[x(t)，u(t)，t]+v(t)

the measurement equation can be specifically expressed as:

$\left(\begin{array}{c}\mathrm{\&theta;}\\ \mathrm{\&phi;}\\ p\\ q\\ r\\ a_x\\ a_y\\ a_z\\ V\end{array}\right)=\left(\begin{array}{c}\mathrm{\&theta;}\\ \mathrm{\&phi;}\\ p\\ q\\ r\\ \frac{1}{m}(T_x+L\sin \mathrm{\&alpha;}-Y\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}-D\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;})\\ \frac{1}{m}(T_y+Y\cos \mathrm{\&beta;}-D\sin \mathrm{\&beta;})\\ \frac{1}{m}(T_z-L\cos \mathrm{\&alpha;}-Y\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-D\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;})\\ V\end{array}\right)+v\left(t\right)$

wherein the noise is measured $v\left(t\right)=[{{\mathrm{\&delta;}}_{\mathrm{\&theta;}},{\mathrm{\&delta;}}_{\mathrm{\&phi;}},{\mathrm{\&delta;}}_p,{\mathrm{\&delta;}}_q,{\mathrm{\&delta;}}_r,{\mathrm{\&delta;}}_{a_x},{\mathrm{\&delta;}}_{a_y},{\mathrm{\&delta;}}_{a_z},{\mathrm{\&delta;}}_V]}^{T^{\&prime;}}$ Caused by the measurement noise of the vacuum velocity calculated by the gyroscope and the linear accelerometer of the inertial device and the GPS;

(c) the state equation and the measurement equation are converted into a continuous linear interference equation:

converting the state equation in the step (a) and the measurement equation in the step (b) into a continuous linear interference equation form by using a Taylor series expansion method, namely f [ X (t), u (t), t)]And h [ X (t), u (t), t)]In-state optimal estimation value

And performing Taylor expansion, and omitting terms with more than two times to obtain a continuous linear interference equation:

$\left\{\begin{array}{c}\mathrm{\&delta;}\stackrel{\&CenterDot;}{X}\left(t\right)=F\left(t\right)\mathrm{\&delta;X}\left(t\right)+G\left(t\right)w\left(t\right)\\ \mathrm{\&delta;Z}\left(t\right)=H\left(t\right)\mathrm{\&delta;X}\left(t\right)+v\left(t\right)\end{array}\right.$

in the formula, $\mathrm{\&delta;X}\left(t\right)=[X\left(t\right)-\hat{X}\left(t\right)],$ $\mathrm{\&delta;Z}\left(t\right)=[Z\left(t\right)-\hat{Z}\left(t\right)],$ $F\left(t\right)=\frac{\&PartialD;f[X\left(t\right),u\left(t\right),t]}{\&PartialD;X\left(t\right)}{|}_{X\left(t\right)=\hat{X}\left(t\right)},$ $H\left(t\right)=\frac{\&PartialD;h[X\left(t\right),u\left(t\right),t]}{\&PartialD;X\left(t\right)}{|}_{X\left(t\right)=\hat{X}\left(t\right)},$ $\hat{Z}\left(t\right)=h[X\left(t\right),u\left(t\right),t]{|}_{X\left(t\right)=\hat{X}\left(t\right)};$

(d) discretizing a continuous linear disturbance equation:

taking the sampling period delta T as T_k+1-t_kDiscretizing the continuous linear disturbance equation obtained in the step (c) to obtain a discrete linear disturbance equation:

$\left\{\begin{array}{c}\mathrm{\&delta;}{X}_{k+1}={\mathrm{\&Phi;}}_{k+1,k}\mathrm{\&delta;}{X}_k+G_k{W}_k\\ \mathrm{\&delta;}{Z}_{k+1}=H_{k+1}\mathrm{\&delta;}{X}_{k+1}+V_{k+1}\end{array}\right.$

in the formula, $\mathrm{\&delta;}{X}_{k+1}=X_{k+1}-{\hat{X}}_{k+1/k},$ $\mathrm{\&delta;}{Z}_{k+1}=Z_{k+1}-{\hat{Z}}_{k+1},$ phi when sampling period DeltaT is small_k+1，k＝I+F(t_k)·ΔT， $G_k=\mathrm{\&Delta;T}\&CenterDot;(I+\frac{F\left(t_k\right)}{2!}\&CenterDot;\mathrm{\&Delta;T})\&CenterDot;G\left(t_k\right),$ $F\left(t_k\right)=\frac{\&PartialD;f[X\left(t_k\right),u\left(t_k\right),t_k]}{\&PartialD;X\left(t_k\right)}{|}_{X\left(t_k\right)={\hat{X}}_k},$ $H_{k+1}=\frac{\&PartialD;h[X\left(t_{k+1}\right),u\left(t_{k+1}\right),t_{k+1}]}{\&PartialD;X\left(t_{k+1}\right)}{|}_{X\left(t_{k+1}\right)={\hat{X}}_{k+1/k}},$

Is the last time instant t_kThe best estimate of the state of the moment in time,for the current time instant, i.e. t_k+1One-step predicted value of the state quantity at the moment;

(e) extended kalman filter equation:

t is obtained according to the step (3)_kAerodynamic force D, Y, L at time, external torque

M, N, engine thrust T_x、T_y、T_zThe coefficient c of the moment equation set obtained in the step (7)₁、c₂、c₃、c₄、c₅、c₆、c₇、c₈、c₉T calculated in step (8)_kTotal mass m of the aircraft at the moment in time, and t_kThree-dimensional wind speed u of a moment_w、v_w、w_wAnd three dimensional rate of change of wind speed

Parameters in combination with t_kOptimal estimation of the state of a time

By using

Finding t_k+1One-step prediction value of time state quantity

t_k+1One-step prediction value of time state quantity

Variance matrix P of_k+1/kPassing through type

Solving is carried out; t is t_k+1Time of day filter gain matrix K_k+1Passing through type $K_{k+1}=P_{k+1/k}{H}_{k+1}^{T^{\&prime;}}[{H_{k+1}{P}_{k+1/k}{H}_{k+1}^{T^{\&prime;}}+R_{k+1}]}^{-1}$ Solving is carried out;

obtaining t according to the step (2)_k+1Three-dimensional ground speed u of time_g、v_g、w_gBy using

To obtain t_k+1Measuring the vacuum speed at the moment, and combining the t obtained in the step (1)_k+1Moment pitch angle theta, roll angle phi, roll angular velocity p, pitch angular velocity q, yaw angular velocity r, and three-dimensional linear acceleration a_x、a_y、a_zParameter composition t_k+1Total quantity of time measurement Z_k+1Combining the solved t_k+1One-step prediction value of time state quantity

And a filter gain matrix K_k+1Using the formula $\mathrm{\&delta;}{\hat{X}}_{k+1}=K_{k+1}\{Z_{k+1}-h[{\hat{X}}_{k+1/k},u\left(t_{k+1}\right),t_{k+1}\left]\right\}$ Is carried out t_k+1Time state quantity improvement value

Solving;

will obtain t_k+1Time state quantity improvement value

And one-step predictionAdd, i.e. toTo obtain t_k+1Optimal estimation of the state of a time

Optimal estimated value

The error variance matrix of (2) can be represented by $P_{k+1}=(I-K_{k+1}{H}_{k+1})P_{k+1/k}{\left({I-K}_{k+1}{H}_{k+1}\right)}^{T^{\&prime;}}+K_{k+1}{R}_{k+1}{K}_{k+1}^{T^{\&prime;}}$ Solving is carried out, thereby realizing t under the flight state with large attack angle_k+1Accurately estimating the moment attack angle and the sideslip angle in real time;

in the extended Kalman filter, initial values of a pitch angle theta, a roll angle phi, a roll angular velocity p, a pitch angular velocity q and a yaw angular velocity r in state quantities are provided by a strapdown inertial navigation system inside a system, and initial values of an attack angle alpha, a sideslip angle beta and a vacuum velocity V in the state quantities are provided by an atmospheric data system outside the system, so that the initial values are formedStatus of state

Initial state estimation error variance matrix P₀Initial value Q of system noise variance array₀Measuring the initial value R of the noise variance matrix₀The initial values of the matrix are directly input from the outside of the system, and the noise matrix Q of the system_kInitial value Q of system noise variance matrix₀Determining and measuring a noise matrix R_k+1By measuring the initial value R of the noise variance matrix₀Determining, in addition, a three-dimensional wind speed u_w、v_w、w_wIs provided by an anemometer external to the system, three-dimensional wind speed rate of change

Are all taken to be zero.

The specific method for calculating the wind speed and the wind speed change rate in the step (10) is as follows:

t obtained according to step (1)_k+1The pitch angle theta, the roll angle phi and the yaw angle psi of the moment are obtained to obtain the body t_k+1Attitude transformation matrix S of time_θψφ：

$S_{\mathrm{\&theta;\&psi;\&phi;}}=\left[\begin{array}{ccc}\cos \mathrm{\&theta;}\cos \mathrm{\&psi;}& \cos \mathrm{\&theta;}\sin \mathrm{\&psi;}& -\sin \mathrm{\&theta;}\\ \sin \mathrm{\&theta;}\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}-\sin \mathrm{\&psi;}\cos \mathrm{\&phi;}& \sin \mathrm{\&theta;}\sin \mathrm{\&psi;}\sin \mathrm{\&phi;}+\cos \mathrm{\&psi;}\cos \mathrm{\&phi;}& \cos \mathrm{\&theta;}\sin \mathrm{\&phi;}\\ \sin \mathrm{\&theta;}\cos \mathrm{\&psi;}\cos \mathrm{\&phi;}+\sin \mathrm{\&psi;}\sin \mathrm{\&phi;}& \sin \mathrm{\&theta;}\sin \mathrm{\&psi;}\cos \mathrm{\&phi;}-\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}& \cos \mathrm{\&theta;}\cos \mathrm{\&phi;}\end{array}\right];$

Combining t obtained in the step (2)_k+1Three-dimensional ground speed u of time_g、v_g、w_gTo obtain t_k+1Three-dimensional ground speed u along x-axis, y-axis and z-axis directions under time machine body coordinate system_e，v_e，w_eComprises the following steps:

$\left[\begin{array}{c}{u}_e\\ v_e\\ w_e\end{array}\right]=S_{\mathrm{\&theta;\&psi;\&phi;}}\left[\begin{array}{c}{u}_g\\ v_g\\ w_g\end{array}\right];$

combining the t estimated in the step (9)_k+1The attack angle alpha, the sideslip angle beta and the vacuum speed V parameter of the moment are converted into a three-dimensional vacuum speed u along the directions of an x axis, a y axis and a z axis under a body coordinate system by using the following formula_b、v_b、w_bI.e. by

$\left[\begin{array}{c}{u}_b\\ v_b\\ w_b\end{array}\right]=\left[\begin{array}{c}V\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}\\ V\sin \mathrm{\&beta;}\\ V\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}\end{array}\right];$

Then t_k+1Three-dimensional wind speed u along x-axis, y-axis and z-axis directions under time machine body coordinate system_w、v_w、w_wThe calculation is made using the following formula:

$\left[\begin{array}{c}{u}_w\\ v_w\\ w_w\end{array}\right]=\left[\begin{array}{c}{u}_e\\ v_e\\ w_e\end{array}\right]-\left[\begin{array}{c}{u}_b\\ v_b\\ w_b\end{array}\right];$

using calculated t_k+1Three-dimensional wind speed u of a moment_w(t_k+1)、v_w(t_k+1)、w_w(t_k+1) Combined with t_kThree-dimensional wind speed u of a moment_w(t_k)、v_w(t_k)、w_w(t_k) By differentiating to find t_k+1Time three-dimensional wind speed rate of change ${\stackrel{\&CenterDot;}{u}}_w\left(t_{k+1}\right),{\stackrel{\&CenterDot;}{v}}_w\left(t_{k+1}\right),{\stackrel{\&CenterDot;}{w}}_w\left(t_{k+1}\right),$ Namely:

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{u}}_w\left(t_{k+1}\right)=\frac{u_w\left(t_{k+1}\right)-u_w\left(t_k\right)}{\mathrm{\&Delta;T}}\\ {\stackrel{\&CenterDot;}{v}}_w\left(t_{k+1}\right)=\frac{v_w\left(t_{k+1}\right)-v_w\left(t_k\right)}{\mathrm{\&Delta;T}}\\ {\stackrel{\&CenterDot;}{w}}_w\left(t_{k+1}\right)=\frac{w_w\left(t_{k+1}\right)-w_w\left(t_k\right)}{\mathrm{\&Delta;T}}\end{array}\right..$

the invention has the advantages and obvious effects that: according to the method, on the premise that the structural appearance of an aircraft is not changed and additional measuring devices and hardware equipment are not added, the output parameters of the existing airborne inertial navigation system, a GPS, a flight control system and a fuel oil measuring system are fully utilized, the angular acceleration, the rotational inertia, the inertia product and the coefficient and the total mass of a moment equation set of the aircraft are solved, an extended Kalman filter is built on the basis of a flight dynamics model, the processes of solving the flight dynamics model and state estimation are integrated, the accurate estimation of the attack angle, the sideslip angle and the vacuum speed at the current moment is realized, the wind speed and the wind speed change rate at the current moment are further solved and input into the extended Kalman filter as feedback information, the real-time accurate estimation of the attack angle and the sideslip angle at the next moment is finished, and the real-time accurate estimation of the attack angle and the sideslip angle is realized in a recursion solving. The method is beneficial to hiding the aircraft, improves the measurement range and measurement precision of the atmospheric data system and the adaptability to complex flight environments, and is relatively economical and easy to implement in engineering. The method has important practical significance for improving the measurement range of the attack angle and the sideslip angle of the new generation of aviation aircrafts in China in the large attack angle flight state and improving the flight performance of the aircrafts.

Drawings

FIG. 1 is a schematic flow diagram of the principle of the method of the present invention;

FIG. 2 is a functional block diagram of a coefficient calculation module of the set of moment equations;

FIG. 3 is a functional block diagram of an extended Kalman filter algorithm module;

FIG. 4 is a functional block diagram of a wind speed and rate of change of wind speed calculation module;

FIG. 5 is a schematic diagram of information interaction between the extended Kalman filter algorithm module and the wind speed and wind speed change rate calculation module;

FIG. 6 is a schematic view of an angle of attack estimation effect curve according to the present invention;

FIG. 7 is a graphical illustration of the side slip angle estimation effect of the present invention;

FIG. 8 is a schematic view of an angle of attack estimation error curve according to the present invention;

FIG. 9 is a schematic diagram of a sideslip angle estimation error curve of the present invention.

Detailed Description

Referring to fig. 1-5, the invention uses the output parameters of the existing airborne inertial navigation system, GPS, flight control system and fuel oil measurement system to calculate the angular acceleration, rotational inertia, inertia product, and coefficient and total mass of the moment equation set of the aircraft, based on the flight dynamics model, by constructing an extended kalman filter, integrating the processes of solving the flight dynamics model and estimating the state, the accurate estimation of the attack angle, the sideslip angle and the vacuum speed at the current moment is realized, and then the wind speed and the wind speed change rate at the current moment are calculated and input into the extended kalman filter as feedback information, the real-time accurate estimation of the attack angle and the sideslip angle at the next moment is completed, and the real-time accurate estimation of the attack angle and the sideslip angle is realized by a recursive solution.

To achieve estimation of the angle of attack and the angle of sideslip in a high angle of attack flight state, the work that needs to be done is:

(1) acquiring data of a strapdown inertial navigation system:

reading three angular velocities, three attitude angles and three linear acceleration information in a strapdown inertial navigation system by using a period delta T, wherein the three angular velocity information are roll angular velocity p, pitch angular velocity q and yaw angular velocity r in the directions of an x axis, a y axis and a z axis under a machine body coordinate system respectively; the three attitude angle information are a pitch angle theta, a roll angle phi and a yaw angle psi respectively; the three pieces of linear acceleration information are respectively the linear acceleration a in the x-axis direction under the body coordinate system_xLinear acceleration a in the y-axis direction_yLinear acceleration a in the z-axis direction_zThe directions of an x axis, a y axis and a z axis of the body coordinate system are forward, rightward and downward respectively;

(2) acquiring GPS system data:

reading the ground speed under the geographic coordinate system output by the GPS system by a period delta T, wherein the directions of the geographic coordinate system are defined as north, east and ground, and the ground speeds in the three directions are u respectively_g、v_g、w_g；

(3) Acquiring flight control system data:

reading the aerodynamic force, the external torque and the engine thrust of the aircraft output by the flight control system at a period delta T, wherein the aerodynamic force of the aircraft is x under an airflow coordinate system_aAxial negative direction, y_aAxial direction and z_aThe components in the negative direction are respectively resistance D, lateral force Y and lift forceL, components of the engine thrust in the directions of the x axis, the y axis and the z axis under a body coordinate system are respectively T_x、T_y、T_zThe components of the external torque in the directions of the x axis, the y axis and the z axis under the body coordinate system are rolling torque respectively

Pitching moment M and yawing moment N;

(4) obtaining the residual fuel quality of the aircraft:

reading the residual fuel mass m output by the fuel measurement system of the aircraft in a period delta T_f；

(5) Angular acceleration calculation:

according to the roll angular velocity p, the pitch angular velocity q and the yaw angular velocity r under the body coordinate system obtained in the step (1), wherein three angular velocity information at the time T are the roll angular velocity p (T), the pitch angular velocity q (T) and the yaw angular velocity r (T), three angular velocity information at the time T + delta T are the roll angular velocity p (T + delta T), the pitch angular velocity q (T + delta T), the yaw angular velocity r (T + delta T), three angular acceleration information in the directions of the x axis, the y axis and the z axis under the body coordinate system at the time T, and the roll angular acceleration

Acceleration of pitch angle

Yaw angular acceleration

Calculated by the differential method, namely:

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{p}\left(t\right)=\frac{p(t+\mathrm{\&Delta;T})-p\left(t\right)}{\mathrm{\&Delta;T}}\\ \stackrel{\&CenterDot;}{q}\left(t\right)=\frac{q(t+\mathrm{\&Delta;T})-q\left(t\right)}{\mathrm{\&Delta;T}}\\ \stackrel{\&CenterDot;}{r}\left(t\right)=\frac{r(t+\mathrm{\&Delta;T})-r\left(t\right)}{\mathrm{\&Delta;T}}\end{array}\right.;$

(6) resolving the rotational inertia and the inertia product of the aircraft:

known aircraft angular motion equation set $\left\{\begin{array}{c}\stackrel{\&CenterDot;}{p}{I}_x-\stackrel{\&CenterDot;}{r}{I}_{\mathrm{xz}}+\mathrm{qr}(I_z-I_y)-\mathrm{pq}{I}_{\mathrm{xz}}=\stackrel{\&OverBar;}{L}\\ \stackrel{\&CenterDot;}{q}{I}_y+\mathrm{pr}(I_x-I_z)+(p^2-r^2)I_{\mathrm{xz}}=M\\ \stackrel{\&CenterDot;}{r}{I}_z-\stackrel{\&CenterDot;}{p}{I}_{\mathrm{xz}}+\mathrm{pq}(I_y-I_x)+\mathrm{qr}{I}_{\mathrm{xz}}=N\end{array}\right.$ Is added to obtain the equation EJ ═ U, where $E=[\stackrel{\&CenterDot;}{p}+\mathrm{pr}-\mathrm{pq},\stackrel{\&CenterDot;}{q}-\mathrm{qr}+\mathrm{pq},\stackrel{\&CenterDot;}{r}+\mathrm{qr}-\mathrm{pr},p^2-r^2-\stackrel{\&CenterDot;}{p}-\stackrel{\&CenterDot;}{r}-\mathrm{pq}+\mathrm{qr}],$ J＝[I_x，I_y，I_z，I_xz]^T′、

T' denotes transpose, I_x、I_y、I_zThe inertia products of the aircraft to the x axis and the y axis are I_xyThe product of inertia about the x-axis and z-axis is I_xzThe product of inertia about the y-axis and the z-axis is I_yzThe aircraft is symmetrical about the plane of symmetry Oxz of the coordinate system of the aircraft body, the product of inertia I_xy＝I_yz＝0，I_xzIs non-zero; from t_aTime begins to t_bIn the time period when the moment is over, continuously acquiring n groups of roll angular velocities p, pitch angular velocities q and yaw angular velocities r by using the step (1), and continuously acquiring n groups of roll torques by using the step (3)

Pitching moment M and yawing moment N, and continuously obtaining N groups of roll angular accelerations in step (5)

Acceleration of pitch angle

Yaw angular accelerationAccording to the obtained n groups of roll angular velocities p and pitchesAngular velocity q, yaw rate r, roll accelerationAcceleration of pitch angle

Yaw angular accelerationRolling torque

The pitching moment M, the yawing moment N and the equation EJ ═ U are calculated by using a least square estimation method to obtain the moment of inertia I of the aircraft_x、I_y、I_zProduct of sum and inertia I_xzWherein n > 4;

(7) and (3) solving the coefficients of the moment equation set:

the moment of inertia I obtained according to the step (6)_x、I_y、I_zProduct of sum and inertia I_xzFor known aircraft moment equations $\left\{\begin{array}{c}\stackrel{\&CenterDot;}{p}=(c_{1}r+c_{2}p)q+c_3\stackrel{\&OverBar;}{L}+c_{4}N\\ \stackrel{\&CenterDot;}{q}={\mathrm{<figure-callout\; id="5"\; label="c"\; filenames="HDA0000121857970000032.png"\; state="\{\{state\}\}">c</figure-callout>}}_5\mathrm{pr}-c_6(p^2-r^2)+c_{7}M\\ \stackrel{\&CenterDot;}{r}=(c_{8}p-c_{2}r)q+c_4\stackrel{\&OverBar;}{L}+c_{9}N\end{array}\right.$ 9 coefficients c in₁、c₂、c₃、c₄、c₅、c₆、c₇、c₈、c₉Carry out a solution in which ${\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="HDA0000121857970000012.png,HDA0000121857970000041.png"\; state="\{\{state\}\}">c</figure-callout>}}_1=\frac{(I_y-I_z)I_z-I_{\mathrm{xz}}^2}{\mathrm{\&Sigma;}},$ $c_2=\frac{(I_x-I_y+I_z)I_{\mathrm{xz}}}{\mathrm{\&Sigma;}},$ $c_3=\frac{I_z}{\mathrm{\&Sigma;}},$ $c_4=\frac{I_{\mathrm{xz}}}{\mathrm{\&Sigma;}},$ ${\mathrm{<figure-callout\; id="5"\; label="c"\; filenames="HDA0000121857970000032.png"\; state="\{\{state\}\}">c</figure-callout>}}_5=\frac{I_z-I_x}{I_y},$ $c_6=\frac{I_{\mathrm{xz}}}{I_y},$ $c_7=\frac{1}{I_y},$ $c_8=\frac{I_x(I_x-I_y)+I_{\mathrm{xz}}^2}{\mathrm{\&Sigma;}},$ $c_9=\frac{I_x}{\mathrm{\&Sigma;}},$ Wherein $\mathrm{\&Sigma;}=I_x{I}_z-I_{\mathrm{xz}}^2;$

(8) The total mass of the aircraft is calculated:

the mass m of the residual fuel oil obtained in the step (4)_fAdding the mass of the aircraft body, the mass of the airborne equipment, the mass of passengers and the mass of weapons to obtain the total mass m of the aircraft

(9) An extended Kalman filter algorithm:

selecting a pitch angle, a roll angle speed, a pitch angle speed, a yaw angle speed, an attack angle, a sideslip angle and a vacuum speed of the aircraft as state quantities according to the flight dynamics model, and further establishing an extended Kalman filter state equation; selecting a pitch angle, a roll angular velocity, a pitch angle velocity, a yaw angular velocity, a three-dimensional linear acceleration and a vacuum velocity as measurement quantities, and further establishing an extended Kalman filter measurement equation; obtaining the current time t according to the step (1) and the step (2)_k+1Measuring information of time, the last time t obtained in step (3)_kAerodynamic force, external torque and engine thrust at the moment, 9 torque equation set coefficients obtained in step (7) and t obtained in step (8)_kTotal mass of the aircraft at time, and t_kThe three-dimensional wind speed and the three-dimensional wind speed change rate at the moment are obtained by using an extended Kalman filter equation_k+1The optimal estimated value of the state quantity at the moment is realized, so that t in the flight state with a large attack angle is realized_k+1Real-time accurate estimation of the angle of attack and sideslip angle at a time, where t_k＞t_b；

(a) Establishment of extended Kalman filter state equation

According to the flight dynamics model, the pitch angle theta, the roll angle phi, the roll angle speed p, the pitch angle speed q, the yaw angle speed r, the attack angle alpha, the sideslip angle beta and the vacuum speed V of the aircraft are selected as state quantities X, namely X is [ theta, phi, p, q, r, alpha, beta, V]^T′And further establishing a state equation of the extended Kalman filter:

$\stackrel{\&CenterDot;}{X}\left(t\right)=f[X\left(t\right),u\left(t\right),t]+G\left(t\right)w\left(t\right)$

the specific form of the equation of state is:

wherein u (t) is control input quantity, and three-dimensional wind speeds along x-axis, y-axis and z-axis directions under a body coordinate system are respectively set as u_w、v_w、w_wThree-dimensional wind speed u_w、v_w、w_wRespectively is

System noise $w\left(t\right)={{[\mathrm{\&delta;}}_{\stackrel{\&OverBar;}{L}},{\mathrm{\&delta;}}_M,{\mathrm{\&delta;}}_N,{\mathrm{\&delta;}}_D,{\mathrm{\&delta;}}_Y,{\mathrm{\&delta;}}_L,{\mathrm{\&delta;}}_{T_x},{\mathrm{\&delta;}}_{T_y},{\mathrm{\&delta;}}_{T_z},{\mathrm{\&delta;}}_{u_w},{\mathrm{\&delta;}}_{v_w},{\mathrm{\&delta;}}_{w_w},{\mathrm{\&delta;}}_{{\stackrel{\&CenterDot;}{u}}_w},{\mathrm{\&delta;}}_{{\stackrel{\&CenterDot;}{v}}_w},{\mathrm{\&delta;}}_{{\stackrel{\&CenterDot;}{w}}_w},{\mathrm{\&delta;}}_m]}^{T^{\&prime;}}$ Respectively by external torque of aircraft

M, N, aerodynamic force D, Y, L, engine thrust T_x、T_y、T_zWind speed u_w、v_w、w_wRate of change of wind speed

And total mass m, g is the acceleration of gravity, and the system noise figure array g (t) is defined as:

$G\left(t\right)=\left[\begin{array}{cccc}{0}_{2\&times;3}& 0_{2\&times;6}& 0_{2\&times;3}& 0_{2\&times;4}\\ A_{3\&times;3}& 0_{3\&times;6}& 0_{3\&times;3}& 0_{3\&times;4}\\ 0_{3\&times;3}& B_{3\&times;6}& C_{3\&times;3}& D_{3\&times;4}\end{array}\right]$

wherein, $A_{3\&times;3}=\left[\begin{array}{ccc}{c}_3& 0& c_4\\ 0& c_7& 0\\ c_4& 0& c_9\end{array}\right]$

$B_{3\&times;6}=\left[\begin{array}{cccccc}0& 0& \frac{-1}{\mathrm{mV}\cos \mathrm{\&beta;}}& \frac{-\sin \mathrm{\&alpha;}}{\mathrm{mV}\cos \mathrm{\&beta;}}& 0& \frac{\cos \mathrm{\&alpha;}}{\mathrm{mV}\cos \mathrm{\&beta;}}\\ 0& \frac{1}{\mathrm{mV}}& 0& \frac{-\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{\mathrm{mV}}& \frac{\cos \mathrm{\&beta;}}{\mathrm{mV}}& \frac{-\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{\mathrm{mV}}\\ \frac{-1}{\mathrm{<figure-callout\; id="0"\; label="m"\; filenames="HDA0000121857970000031.png,HDA0000121857970000032.png"\; state="\{\{state\}\}">m</figure-callout>}}& 0& 0& \frac{\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}}{m}& \frac{\sin \mathrm{\&beta;}}{m}& \frac{\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}}{m}\end{array}\right]$

$C_{3\&times;3}=\left[\begin{array}{ccc}\frac{q\cos \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}& \frac{-p\cos \mathrm{\&alpha;}-r\sin \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}& \frac{q\sin \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}\\ \frac{-r\cos \mathrm{\&beta;}-q\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}& \frac{-r\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}+p\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}& \frac{p\cos \mathrm{\&beta;}+q\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}\\ -r\sin \mathrm{\&beta;}+q\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}& r\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}-p\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}& p\sin \mathrm{\&beta;}-q\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}\end{array}\right]$

$D_{3\&times;4}=\left[\begin{array}{cccc}\frac{\sin \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}& 0& \frac{-\cos \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}& \frac{L+T_x\sin \mathrm{\&alpha;}-T_z\cos \mathrm{\&alpha;}}{m^{2}V\cos \mathrm{\&beta;}}\\ \frac{\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}& \frac{-\cos \mathrm{\&beta;}}{V}& \frac{\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}& \frac{-Y+T_x\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}-T_y\cos \mathrm{\&beta;}+T_z\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{m^{2}V}\\ -\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}& -\sin \mathrm{\&beta;}& -\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}& \frac{D-T_x\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}-T_y\sin \mathrm{\&beta;}-T_z\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}}{m^2}\end{array}\right]$

(b) establishment of extended Kalman filter measurement equation

Three-dimensional ground speed u obtained from GPS_g、v_g、w_gBy using

The calculated V is used as the vacuum speed measurement and is combined with the pitch angle theta, the roll angle phi, the roll angle speed p, the pitch angle speed q, the yaw angle speed r and the three-dimensional linear acceleration a which are obtained by the strapdown inertial navigation system_x、a_y、a_zThe quantity of the total composition Z, i.e. Z ═ θ, φ, p, q, r, a_x，a_y，a_z，V]^T′And further establishing an extended Kalman filter measurement equation:

Z(t)＝h[X(t)，u(t)，t]+v(t)

the measurement equation can be specifically expressed as:

$\left(\begin{array}{c}\mathrm{\&theta;}\\ \mathrm{\&phi;}\\ p\\ q\\ r\\ a_x\\ a_y\\ a_z\\ V\end{array}\right)=\left(\begin{array}{c}\mathrm{\&theta;}\\ \mathrm{\&phi;}\\ p\\ q\\ r\\ \frac{1}{m}(T_x+L\sin \mathrm{\&alpha;}-Y\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}-D\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;})\\ \frac{1}{m}(T_y+Y\cos \mathrm{\&beta;}-D\sin \mathrm{\&beta;})\\ \frac{1}{m}(T_z-L\cos \mathrm{\&alpha;}-Y\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-D\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;})\\ V\end{array}\right)+v\left(t\right)$

wherein the noise is measured $v\left(t\right)=[{{\mathrm{\&delta;}}_{\mathrm{\&theta;}},{\mathrm{\&delta;}}_{\mathrm{\&phi;}},{\mathrm{\&delta;}}_p,{\mathrm{\&delta;}}_q,{\mathrm{\&delta;}}_r,{\mathrm{\&delta;}}_{a_x},{\mathrm{\&delta;}}_{a_y},{\mathrm{\&delta;}}_{a_z},{\mathrm{\&delta;}}_V]}^{T^{\&prime;}}$ Caused by the measurement noise of the vacuum velocity calculated by the gyroscope and the linear accelerometer of the inertial device and the GPS;

(c) the state equation and the measurement equation are converted into a continuous linear interference equation:

converting the state equation in the step (a) and the measurement equation in the step (b) into a continuous linear interference equation form by using a Taylor series expansion method, namely f [ X (t), u (t), t)]And h [ X (t), u (t), t)]In-state optimal estimation value

And performing Taylor expansion, and omitting terms with more than two times to obtain a continuous linear interference equation:

$\left\{\begin{array}{c}\mathrm{\&delta;}\stackrel{\&CenterDot;}{X}\left(t\right)=F\left(t\right)\mathrm{\&delta;X}\left(t\right)+G\left(t\right)w\left(t\right)\\ \mathrm{\&delta;Z}\left(t\right)=H\left(t\right)\mathrm{\&delta;X}\left(t\right)+v\left(t\right)\end{array}\right.$

in the formula, $\mathrm{\&delta;X}\left(t\right)=[X\left(t\right)-\hat{X}\left(t\right)],$ $\mathrm{\&delta;Z}\left(t\right)=[Z\left(t\right)-\hat{Z}\left(t\right)],$ $F\left(t\right)=\frac{\&PartialD;f[X\left(t\right),u\left(t\right),t]}{\&PartialD;X\left(t\right)}{|}_{X\left(t\right)=\hat{X}\left(t\right)},$ $H\left(t\right)=\frac{\&PartialD;h[X\left(t\right),u\left(t\right),t]}{\&PartialD;X\left(t\right)}{|}_{X\left(t\right)=\hat{X}\left(t\right)},$ $\hat{Z}\left(t\right)=h[X\left(t\right),u\left(t\right),t]{|}_{X\left(t\right)=\hat{X}\left(t\right)};$

(d) discretizing a continuous linear disturbance equation:

taking the sampling period delta T as T_k+1-t_kDiscretizing the continuous linear disturbance equation obtained in the step (c) to obtain a discrete linear disturbance equation:

$\left\{\begin{array}{c}\mathrm{\&delta;}{X}_{k+1}={\mathrm{\&Phi;}}_{k+1,k}\mathrm{\&delta;}{X}_k+G_k{W}_k\\ \mathrm{\&delta;}{Z}_{k+1}=H_{k+1}\mathrm{\&delta;}{X}_{k+1}+V_{k+1}\end{array}\right.$

in the formula, $\mathrm{\&delta;}{X}_{k+1}=X_{k+1}-{\hat{X}}_{k+1/k},$ $\mathrm{\&delta;}{Z}_{k+1}=Z_{k+1}-{\hat{Z}}_{k+1},$ phi when sampling period DeltaT is small_k+1，k＝I+F(t_k)·ΔT， $G_k=\mathrm{\&Delta;T}\&CenterDot;(I+\frac{F\left(t_k\right)}{2!}\&CenterDot;\mathrm{\&Delta;T})\&CenterDot;G\left(t_k\right),$ $F\left(t_k\right)=\frac{\&PartialD;f[X\left(t_k\right),u\left(t_k\right),t_k]}{\&PartialD;X\left(t_k\right)}{|}_{X\left(t_k\right)={\hat{X}}_k},$ $H_{k+1}=\frac{\&PartialD;h[X\left(t_{k+1}\right),u\left(t_{k+1}\right),t_{k+1}]}{\&PartialD;X\left(t_{k+1}\right)}{|}_{X\left(t_{k+1}\right)={\hat{X}}_{k+1/k}},$

Is the last time instant t_kThe best estimate of the state of the moment in time,

for the current time instant, i.e. t_k+1One-step predicted value of the state quantity at the moment;

(e) extended kalman filter equation:

t is obtained according to the step (3)_kAerodynamic force D, Y, L at time, external torque

M, N, engine thrust T_x、T_y、T_zThe coefficient c of the moment equation set obtained in the step (7)₁、c₂、c₃、c₄、c₅、c₆、c₇、c₈、c₉T calculated in step (8)_kTotal mass m of the aircraft at the moment in time, and t_kThree-dimensional wind speed u of a moment_w、v_w、w_wAnd three dimensional rate of change of wind speed

Parameters in combination with t_kOptimal estimation of the state of a time

By usingFinding t_k+1One-step prediction value of time state quantity

t_k+1One-step prediction value of time state quantity

Variance matrix P of_k+1/kPassing through type

Solving is carried out; t is t_k+1Time of day filter gain matrix K_k+1Passing through type $K_{k+1}=P_{k+1/k}{H}_{k+1}^{T^{\&prime;}}[{H_{k+1}{P}_{k+1/k}{H}_{k+1}^{T^{\&prime;}}+R_{k+1}]}^{-1}$ Solving is carried out;

obtaining t according to the step (2)_k+1Three-dimensional ground speed u of time_g、v_g、w_gBy using

To obtain t_k+1Measuring the vacuum speed at the moment, and combining the t obtained in the step (1)_k+1Moment pitch angle theta, roll angle phi, roll angular velocity p, pitch angular velocity q, yaw angular velocity r, and three-dimensional linear acceleration a_x、a_y、a_zParameter composition t_k+1Total quantity of time measurement Z_k+1Combining the solved t_k+1One-step prediction value of time state quantity

And a filter gain matrix K_k+1Using the formula $\mathrm{\&delta;}{\hat{X}}_{k+1}=K_{k+1}\{Z_{k+1}-h[{\hat{X}}_{k+1/k},u\left(t_{k+1}\right),t_{k+1}\left]\right\}$ Is carried out t_k+1Time state quantity improvement valueSolving;

will obtain t_k+1Time state quantity improvement value

And one-step prediction

Add, i.e. to

To obtain t_k+1Optimal estimation of the state of a time

Optimal estimated value

The error variance matrix of (2) can be represented by $P_{k+1}=(I-K_{k+1}{H}_{k+1})P_{k+1/k}{\left({I-K}_{k+1}{H}_{k+1}\right)}^{T^{\&prime;}}+K_{k+1}{R}_{k+1}{K}_{k+1}^{T^{\&prime;}}$ Solving is carried out, thereby realizing t under the flight state with large attack angle_k+1Accurately estimating the moment attack angle and the sideslip angle in real time;

in the extended Kalman filter, initial values of a pitch angle theta, a roll angle phi, a roll angular velocity p, a pitch angular velocity q and a yaw angular velocity r in state quantities are provided by a strapdown inertial navigation system in the system, and an attack angle alpha, a sideslip angle beta and a true yaw angle r in the state quantitiesThe initial value of airspeed V is provided by an atmospheric data system external to the system, constituting an initial state

Initial state estimation error variance matrix P₀Initial value Q of system noise variance array₀Measuring the initial value R of the noise variance matrix₀The initial values of the matrix are directly input from the outside of the system, and the noise matrix Q of the system_kInitial value Q of system noise variance matrix₀Determining and measuring a noise matrix R_k+1By measuring the initial value R of the noise variance matrix₀Determining, in addition, a three-dimensional wind speed u_w、v_w、w_wIs provided by an anemometer external to the system, three-dimensional wind speed rate of change

Are all taken to be zero.

(10) Resolving the wind speed and the wind speed change rate:

obtaining t according to step (1)_k+1Pitch angle, roll angle, yaw angle at the moment and t obtained in step (2)_k+1Combining the three-dimensional ground speed parameter of the moment with the t estimated in the step (9)_k+1The attack angle, the sideslip angle and the vacuum speed at the moment are calculated to obtain t_k+1Three dimensional wind speed at time, combined with t_kThe three-dimensional wind speed at the moment is calculated to obtain t_k+1The three-dimensional wind speed rate of change at the moment;

t obtained according to step (1)_k+1The pitch angle theta, the roll angle phi and the yaw angle psi of the moment are obtained to obtain the body t_k+1Attitude transformation matrix S of time_θψφ：

$S_{\mathrm{\&theta;\&psi;\&phi;}}=\left[\begin{array}{ccc}\cos \mathrm{\&theta;}\cos \mathrm{\&psi;}& \cos \mathrm{\&theta;}\sin \mathrm{\&psi;}& -\sin \mathrm{\&theta;}\\ \sin \mathrm{\&theta;}\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}-\sin \mathrm{\&psi;}\cos \mathrm{\&phi;}& \sin \mathrm{\&theta;}\sin \mathrm{\&psi;}\sin \mathrm{\&phi;}+\cos \mathrm{\&psi;}\cos \mathrm{\&phi;}& \cos \mathrm{\&theta;}\sin \mathrm{\&phi;}\\ \sin \mathrm{\&theta;}\cos \mathrm{\&psi;}\cos \mathrm{\&phi;}+\sin \mathrm{\&psi;}\sin \mathrm{\&phi;}& \sin \mathrm{\&theta;}\sin \mathrm{\&psi;}\cos \mathrm{\&phi;}-\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}& \cos \mathrm{\&theta;}\cos \mathrm{\&phi;}\end{array}\right];$

Combining t obtained in the step (2)_k+1Three-dimensional ground speed u of time_g、v_g、w_gTo obtain t_k+1Three-dimensional ground speed u along x-axis, y-axis and z-axis directions under time machine body coordinate system_e，v_e，w_eComprises the following steps:

$\left[\begin{array}{c}{u}_e\\ v_e\\ w_e\end{array}\right]=S_{\mathrm{\&theta;\&psi;\&phi;}}\left[\begin{array}{c}{u}_g\\ v_g\\ w_g\end{array}\right];$

bonding ofT estimated in step (9)_k+1The attack angle alpha, the sideslip angle beta and the vacuum speed V parameter of the moment are converted into a three-dimensional vacuum speed u along the directions of an x axis, a y axis and a z axis under a body coordinate system by using the following formula_b、v_b、w_bI.e. by

$\left[\begin{array}{c}{u}_b\\ v_b\\ w_b\end{array}\right]=\left[\begin{array}{c}V\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}\\ V\sin \mathrm{\&beta;}\\ V\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}\end{array}\right];$

Then t_k+1Three-dimensional wind speed u along x-axis, y-axis and z-axis directions under time machine body coordinate system_w、v_w、w_wThe calculation is made using the following formula:

$\left[\begin{array}{c}{u}_w\\ v_w\\ w_w\end{array}\right]=\left[\begin{array}{c}{u}_e\\ v_e\\ w_e\end{array}\right]-\left[\begin{array}{c}{u}_b\\ v_b\\ w_b\end{array}\right];$

using calculated t_k+1Three-dimensional wind speed u of a moment_w(t_k+1)、v_w(t_k+1)、w_w(t_k+1) Combined with t_kThree-dimensional wind speed u of a moment_w(t_k)、v_w(t_k)、w_w(t_k) By differentiating to find t_k+1Time three-dimensional wind speed rate of change ${\stackrel{\&CenterDot;}{u}}_w\left(t_{k+1}\right),{\stackrel{\&CenterDot;}{v}}_w\left(t_{k+1}\right),{\stackrel{\&CenterDot;}{w}}_w\left(t_{k+1}\right),$ Namely:

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{u}}_w\left(t_{k+1}\right)=\frac{u_w\left(t_{k+1}\right)-u_w\left(t_k\right)}{\mathrm{\&Delta;T}}\\ {\stackrel{\&CenterDot;}{v}}_w\left(t_{k+1}\right)=\frac{v_w\left(t_{k+1}\right)-v_w\left(t_k\right)}{\mathrm{\&Delta;T}}\\ {\stackrel{\&CenterDot;}{w}}_w\left(t_{k+1}\right)=\frac{w_w\left(t_{k+1}\right)-w_w\left(t_k\right)}{\mathrm{\&Delta;T}}\end{array}\right..$

(11) wind speed and wind speed change rate feedback step:

t obtained in the step (10)_k+1Feeding back the three-dimensional wind speed and the three-dimensional wind speed change rate at the moment to the extended Kalman filtering module algorithm module for completing the next moment, namely t in the step (9)_k+2At the moment of attack and sideslipAnd (6) estimating.

Therefore, the real-time accurate estimation of the attack angle and the sideslip angle can be realized through a recursion solving mode.

FIG. 1 is a basic principle block diagram of the method of the present invention, and the attack angle/sideslip angle calculation system finally realizes the accurate estimation of the attack angle and the sideslip angle by using parameters provided by a strapdown inertial navigation system, a GPS system, a flight control system and a fuel oil measurement system through a moment equation set coefficient calculation module, a wind speed and wind speed change rate calculation module, a total mass calculation module and an extended Kalman filtering algorithm module, wherein solid arrows between the modules represent a basic logical connection relationship.

Fig. 2 is a further refinement of the coefficient calculation module of the moment equation set in fig. 1, which represents a specific flow of coefficient calculation of the moment equation set, and more clearly shows the source and transmission process of information.

Fig. 3 is a further refinement of the extended kalman filter algorithm block in fig. 1, which represents a specific process of calculating the attack angle, the sideslip angle and the vacuum velocity, and more clearly shows the source of information and the main components of the extended kalman filter.

FIG. 4 is a further refinement of the wind speed and wind speed change rate calculation module shown in FIG. 1, which more clearly shows the calculation process and information source of the wind speed and wind speed change rate at the current time.

FIG. 5 is a further refinement of the information interaction between the extended Kalman filter algorithm module and the wind speed and wind speed change rate calculation module of FIG. 1, more clearly illustrating the process of mutual transmission and utilization of information between the two modules at different times, wherein

Is an initial value of a state quantity, t₁In order to start the extended kalman filter,

is t₁Estimate of the state quantity at time, t₂Is t₁The next moment in time of the first time,

is t₂An estimate of the state quantity at the time of day,

is t_k+1An estimate of the state quantity at the time of day,

is t_k+2An estimate of the state quantity at the time.

FIG. 6 is a curve of the effect of the estimation of the angle of attack of the method of the present invention, the horizontal axis represents time, the vertical axis represents the angle of attack, the variation range of the flight angle of attack is-41 to-75 degrees, the angle of attack estimated by the present invention is basically coincident with the curve of the true angle of attack, the goodness of fit is good, the present invention can be applied to not only realize the accurate estimation of the angle of attack under the conventional flight angle of attack, but also have good tracking effect for the flight with large angle of attack.

FIG. 7 is a curve of the effect of estimating the sideslip angle by the method of the present invention, wherein the horizontal axis represents time, the vertical axis represents the sideslip angle, the variation range of the sideslip angle is-12 to 22 degrees, the sideslip angle estimated by the method of the present invention is basically coincident with the curve of the real sideslip angle, and the sideslip angle under the flight state with a large attack angle can be accurately estimated by applying the method of the present invention.

FIG. 8 is an estimation error curve of the angle of attack of the method of the present invention, the horizontal axis represents time, the vertical axis represents estimation error of the angle of attack, the mean value of the estimated error of the angle of attack using the method of the present invention is 0.31 degrees, the standard deviation of the estimated error of the angle of attack is 0.20 degrees, and the maximum error of the estimated error of the angle of attack is controlled within 1.2 degrees.

FIG. 9 is a curve of the estimation error of the sideslip angle of the method of the present invention, the horizontal axis represents time, the vertical axis represents the estimation error of the sideslip angle, the mean value of the estimated sideslip angle error by applying the present invention is 0.05 degrees, the standard deviation of the estimated sideslip angle error is 0.33 degrees, and the estimation error of the sideslip angle is controlled within 1.2 degrees.

DETAILED DESCRIPTION OF EMBODIMENT (S) OF INVENTION

In the embodiment, based on flight data simulated by flight simulation software X-Plane, the method for estimating the attack angle and the sideslip angle in the flight state with the large attack angle is tested, and a beneficial conclusion is obtained.

Simulation test effect based on X-Plane system: as can be seen from the attack angle and sideslip angle estimation effect curve diagrams 6 and 7 and the estimation error curve diagrams 8 and 9, the method not only can realize accurate estimation of the attack angle and the sideslip angle in conventional flight, but also can have better tracking effect on large attack angle flight, wherein the average values of the estimation errors of the attack angle and the sideslip angle are respectively 0.31 degrees and 0.05 degrees, the standard deviations of the estimation errors are respectively 0.20 degrees and 0.33 degrees, and although the errors of the attack angle and the sideslip angle estimated in the large attack angle flight phase are slightly increased, the errors of the estimated attack angle and the sideslip angle can be controlled within 1.2 degrees. The invention has strong engineering application value.

Claims (4)

1. An atmospheric attack angle and sideslip angle estimation method in a large attack angle flight state is characterized by comprising the following steps: parameters output using equipment onboard an aircraft, including:

firstly, outputting an attitude angle, an angular velocity and a linear acceleration by a strapdown inertial navigation system;

the ground speed output by the GPS system;

aerodynamic force, external torque and engine thrust output by the flight control system;

fourthly, the residual fuel quality output by the fuel measuring system;

solving angular acceleration through a differential method according to the angular velocity;

solving the rotational inertia and the inertia product of the aircraft by the angular velocity, the angular acceleration and the external moment through a least square estimation method, and further obtaining the coefficient of a moment equation set;

adding the mass of the residual fuel oil with the known mass of an aircraft body, the mass of airborne equipment, the mass of passengers and the mass of weapons to obtain the total mass of the aircraft;

constructing an extended Kalman filter according to the flight dynamics model, and realizing accurate estimation of an attack angle, a sideslip angle and a vacuum speed at the current moment in a large attack angle flight state by using an extended Kalman filtering algorithm;

calculating the wind speed at the current moment according to the attitude angle and the ground speed at the current moment and the estimated attack angle, the sideslip angle and the vacuum speed at the current moment, and calculating the wind speed change rate at the current moment by combining the wind speed at the previous moment;

and feeding the obtained wind speed and the wind speed change rate at the current moment back to the extended Kalman filter mode algorithm module for finishing the real-time accurate estimation of the attack angle and the sideslip angle at the next moment, and realizing the real-time accurate estimation of the attack angle and the sideslip angle in a recursion solving mode.

2. The method for estimating the atmospheric angle of attack and the sideslip angle under the flight condition with the large angle of attack according to claim 1, comprising the steps of:

(1) reading three angular velocities, three attitude angles and three linear acceleration information in a strapdown inertial navigation system by using a period delta T, wherein the three angular velocity information are roll angular velocity p, pitch angular velocity q and yaw angular velocity r in the directions of an x axis, a y axis and a z axis under a machine body coordinate system respectively; the three attitude angle information are a pitch angle theta, a roll angle phi and a yaw angle psi respectively; the three pieces of linear acceleration information are respectively the linear acceleration a in the x-axis direction under the body coordinate system_xLinear acceleration a in the y-axis direction_yLinear acceleration a in the z-axis direction_zThe directions of an x axis, a y axis and a z axis of the body coordinate system are forward, rightward and downward respectively;

(2) reading the ground speed under the geographic coordinate system output by the GPS system by a period delta T, wherein the directions of the geographic coordinate system are defined as north, east and ground, and the ground speeds in the three directions are u respectively_g、v_g、w_g；

(3) Reading the aerodynamic force, the external torque and the engine thrust of the aircraft output by the flight control system at a period delta T, wherein the aerodynamic force of the aircraft is x under an airflow coordinate system_aAxial negative direction, y_aAxial direction and z_aThe components in the negative direction of the shaft are respectively resistance D, lateral force Y and lift L, and the components in the directions of the x-axis, the Y-axis and the z-axis of the thrust of the engine under a body coordinate system are respectively T_x、T_y、T_zThe components of the external torque in the directions of the x axis, the y axis and the z axis under the body coordinate system are rolling torque respectively

Pitching moment M and yawing moment N;

(4) reading the residual fuel mass m output by the fuel measurement system of the aircraft in a period delta T_f；

(5) According to the roll angular velocity p, the pitch angular velocity q and the yaw angular velocity r under the body coordinate system obtained in the step (1), wherein three angular velocity information at the time T are the roll angular velocity p (T), the pitch angular velocity q (T) and the yaw angular velocity r (T), and three angular velocity information at the time T + delta T are the roll angular velocity p (T + delta T), the pitch angular velocity q (T + delta T), the yaw angular velocity r (T + delta T), and three angular acceleration information in the directions of the x axis, the y axis and the z axis under the body coordinate system at the time T, namely the roll angular acceleration

Acceleration of pitch angle

Yaw angular accelerationCalculated by the differential method, namely:

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{p}\left(t\right)=\frac{p(t+\mathrm{\&Delta;T})-p\left(t\right)}{\mathrm{\&Delta;T}}\\ \stackrel{\&CenterDot;}{q}\left(t\right)=\frac{q(t+\mathrm{\&Delta;T})-q\left(t\right)}{\mathrm{\&Delta;T}}\\ \stackrel{\&CenterDot;}{r}\left(t\right)=\frac{r(t+\mathrm{\&Delta;T})-r\left(t\right)}{\mathrm{\&Delta;T}}\end{array}\right.;$

(6) known aircraft angular motion equation set $\left\{\begin{array}{c}\stackrel{\&CenterDot;}{p}{I}_x-\stackrel{\&CenterDot;}{r}{I}_{\mathrm{xz}}+\mathrm{qr}(I_z-I_y)-\mathrm{pq}{I}_{\mathrm{xz}}=\stackrel{\&OverBar;}{L}\\ \stackrel{\&CenterDot;}{q}{I}_y+\mathrm{pr}(I_x-I_z)+(p^2-r^2)I_{\mathrm{xz}}=M\\ \stackrel{\&CenterDot;}{r}{I}_z-\stackrel{\&CenterDot;}{p}{I}_{\mathrm{xz}}+\mathrm{pq}(I_y-I_x)+\mathrm{qr}{I}_{\mathrm{xz}}=N\end{array}\right.$ Is added to give the equation EJ = U, where $E=[\stackrel{\&CenterDot;}{p}+\mathrm{pr}-\mathrm{pq},\stackrel{\&CenterDot;}{q}-\mathrm{qr}+\mathrm{pq},\stackrel{\&CenterDot;}{r}+\mathrm{qr}-\mathrm{pr},p^2-r^2-\stackrel{\&CenterDot;}{p}-\stackrel{\&CenterDot;}{r}-\mathrm{pq}+\mathrm{qr}],$ J=[I_x,I_y,I_z,I_xz]^T'、

T' denotes transpose, I_x、I_y、I_zThe inertia products of the aircraft to the x axis and the y axis are I_xyThe product of inertia about the x-axis and z-axis is I_xzThe product of inertia about the y-axis and the z-axis is I_yzThe aircraft is symmetrical about the plane of symmetry Oxz of the coordinate system of the aircraft body, the product of inertia I_xy=I_yz=0，I_xzIs non-zero; from t_aTime begins to t_bIn the time period when the moment is over, continuously acquiring n groups of roll angular velocities p, pitch angular velocities q and yaw angular velocities r by using the step (1), and continuously acquiring n groups of roll torques by using the step (3)

Pitching moment M and yawing moment N, and continuously obtaining N groups of roll angular accelerations in step (5)

Acceleration of pitch angle

Yaw angular acceleration

Obtaining n groups of roll angular velocity p, pitch angular velocity q, yaw angular velocity r and roll angular accelerationAcceleration of pitch angle

Yaw angular accelerationRolling torque

Pitching moment M, yawing moment N and equation EJ = U, and solving the rotational inertia I of the aircraft by using a least square estimation method_x、I_y、I_zProduct of sum and inertia I_xzWherein n is>4；

(7) The moment of inertia I obtained according to the step (6)_x、I_y、I_zProduct of sum and inertia I_xzFor known aircraft moment equations $\left\{\begin{array}{c}\stackrel{\&CenterDot;}{p}=(c_{1}r+c_{2}p)q+c_3\stackrel{\&OverBar;}{L}+c_{4}N\\ \stackrel{\&CenterDot;}{q}=c_5\mathrm{pr}-c_6(p^2-r^2)+c_{7}M\\ \stackrel{\&CenterDot;}{r}=(c_{8}p-c_{2}r)q+c_4\stackrel{\&OverBar;}{L}+c_{9}N\end{array}\right.$ 9 coefficients c in₁、c₂、c₃、c₄、c₅、c₆、c₇、c₈、c₉Carry out a solution in which $c_1=\frac{(I_y-I_z)I_z-I_{\mathrm{xz}}^2}{\mathrm{\&Sigma;}},$ $c_2=\frac{(I_x-I_y+I_z)I_{\mathrm{xz}}}{\mathrm{\&Sigma;}},$ $c_3=\frac{I_z}{\mathrm{\&Sigma;}},$ $c_4=\frac{I_{\mathrm{xz}}}{\mathrm{\&Sigma;}},$ $c_5=\frac{I_z-I_x}{I_y},$ $c_6=\frac{I_{\mathrm{xz}}}{I_y},$ $c_7=\frac{1}{I_y},$ $c_8=\frac{I_x(I_x-I_y)+I_{\mathrm{xz}}^2}{\mathrm{\&Sigma;}},$ $c_9=\frac{I_x}{\mathrm{\&Sigma;}},$ Wherein $\mathrm{\&Sigma;}=I_x{I}_z-I_{\mathrm{xz}}^2;$

(8) The mass m of the residual fuel oil obtained in the step (4)_fAdding the mass of the aircraft body, the mass of airborne equipment, the mass of passengers and the mass of weapons to obtain the total mass m of the aircraft;

(9) selecting a pitch angle, a roll angle speed, a pitch angle speed, a yaw angle speed, an attack angle, a sideslip angle and a vacuum speed of the aircraft as state quantities according to the flight dynamics model, and further establishing an extended Kalman filter state equation; selecting a pitch angle, a roll angular velocity, a pitch angle velocity, a yaw angular velocity, a three-dimensional linear acceleration and a vacuum velocity as measurement quantities, and further establishing an extended Kalman filter measurement equation; obtaining the current time t according to the step (1) and the step (2)_k+1Measuring information of time, the last time t obtained in step (3)_kAerodynamic force, external torque and engine thrust at the moment, 9 torque equation set coefficients obtained in step (7) and t obtained in step (8)_kTotal mass of the aircraft at time, and t_kThe three-dimensional wind speed and the three-dimensional wind speed change rate at the moment are obtained by using an extended Kalman filter equation_k+1The optimal estimated value of the state quantity at the moment is realized, so that t in the flight state with a large attack angle is realized_k+1Real-time accurate estimation of the angle of attack and sideslip angle at a time, where t_k>t_b；

(10) Obtaining t according to step (1)_k+1Pitch angle, roll angle, yaw angle at the moment and t obtained in step (2)_k+1Combining the three-dimensional ground speed parameter of the moment with the t estimated in the step (9)_k+1The attack angle, the sideslip angle and the vacuum speed at the moment are calculated to obtain t_k+1Three dimensional wind speed at time, combined with t_kThe three-dimensional wind speed at the moment is calculated to obtain t_k+1The three-dimensional wind speed rate of change at the moment;

(11) t obtained in the step (10)_k+1Feeding back the three-dimensional wind speed and the three-dimensional wind speed change rate at the moment to the extended Kalman filtering module algorithm module for completing the next moment, namely t in the step (9)_k+2And estimating the moment attack angle and the sideslip angle.

3. The method according to claim 2, wherein the method for estimating the atmospheric angle of attack and the sideslip angle in the flight state with the large angle of attack comprises: the extended kalman filter algorithm in the step (9) specifically comprises the following steps:

(a) establishment of extended Kalman filter state equation

Selecting a pitch angle theta, a roll angle phi, a roll angle speed p, a pitch angle speed q, a yaw angle speed r, an attack angle alpha, a sideslip angle beta and a vacuum speed V of the aircraft as state quantities X according to a flight dynamics model, wherein X = [ theta, phi, p, q, r, alpha, beta, V ] is]^T' further establishing a state equation of the extended Kalman filter:

$\stackrel{\&CenterDot;}{X}\left(t\right)=f[X\left(t\right),u\left(t\right),t]+G\left(t\right)w\left(t\right)$

the specific form of the equation of state is:

$\left(\begin{array}{c}\stackrel{\&CenterDot;}{\mathrm{\&theta;}}\\ \stackrel{\&CenterDot;}{\mathrm{\&phi;}}\\ \stackrel{\&CenterDot;}{p}\\ \stackrel{\&CenterDot;}{q}\\ \stackrel{\&CenterDot;}{r}\\ \stackrel{\&CenterDot;}{\mathrm{\&alpha;}}\\ \stackrel{\&CenterDot;}{\mathrm{\&beta;}}\\ \stackrel{\&CenterDot;}{V}\end{array}\right)=\left(\begin{array}{c}q\cos \mathrm{\&phi;}-r\sin \mathrm{\&phi;}\\ p+(r\cos \mathrm{\&phi;}+q\sin \mathrm{\&phi;})\tan \mathrm{\&theta;}\\ (c_{1}r+c_{2}p)q+c_3\stackrel{\&OverBar;}{L}+c_{4}N\\ c_5\mathrm{pr}-c_6(p^2-r^2)+c_{7}M\\ (c_{8}p-c_{2}r)q+c_4\stackrel{\&OverBar;}{L}+c_{9}N\\ \frac{1}{\mathrm{mV}\cos \mathrm{\&beta;}}[-T_x\sin \mathrm{\&alpha;}+T_z\cos \mathrm{\&alpha;}-L+\mathrm{mV}(-p\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}+q\cos \mathrm{\&beta;}-r\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;})+\\ \mathrm{mg}(\sin \mathrm{\&alpha;}\sin \mathrm{\&theta;}+\cos \mathrm{\&alpha;}\cos \mathrm{\&phi;}\cos \mathrm{\&theta;})+m({\mathrm{qw}}_w-{\mathrm{rv}}_w+{\stackrel{\&CenterDot;}{u}}_w)\sin \mathrm{\&alpha;}-m({\mathrm{pv}}_w-{\mathrm{qu}}_w+{\stackrel{\&CenterDot;}{w}}_w)\cos \mathrm{\&alpha;}]\\ \frac{1}{\mathrm{mV}}[-T_x\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}+T_y\cos \mathrm{\&beta;}-T_z\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}+Y-\mathrm{mV}(-p\sin \mathrm{\&alpha;}+r\cos \mathrm{\&alpha;})+\\ \mathrm{mg}(\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}\sin \mathrm{\&theta;}+\cos \mathrm{\&beta;}\sin \mathrm{\&phi;}\cos \mathrm{\&theta;}-\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}\cos \mathrm{\&phi;}\cos \mathrm{\&theta;})+\\ m({\mathrm{qw}}_w-{\mathrm{rv}}_w+{\stackrel{\&CenterDot;}{u}}_w)\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}+m({\mathrm{pw}}_w-{\mathrm{ru}}_w-{\stackrel{\&CenterDot;}{v}}_w)\cos \mathrm{\&beta;}+m({\mathrm{pv}}_w-{\mathrm{qu}}_w+{\stackrel{\&CenterDot;}{w}}_w)\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}]\\ \frac{1}{m}[T_x\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}+T_y\sin \mathrm{\&beta;}+T_z\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}-D+\mathrm{mg}(-\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}\sin \mathrm{\&theta;}+\\ \sin \mathrm{\&beta;}\sin \mathrm{\&phi;}\cos \mathrm{\&theta;}+\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}\cos \mathrm{\&phi;}\cos \mathrm{\&theta;})-m({\mathrm{qw}}_w-{\mathrm{rv}}_w+{\stackrel{\&CenterDot;}{u}}_w)\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}+\\ m({\mathrm{pw}}_w-{\mathrm{ru}}_w-{\stackrel{\&CenterDot;}{v}}_w)\sin \mathrm{\&beta;}-m({\mathrm{pv}}_w-{\mathrm{qu}}_w+{\stackrel{\&CenterDot;}{w}}_w)\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}]\end{array}\right)+G\left(t\right)w\left(t\right)$

wherein,u (t) is control input quantity, and three-dimensional wind speeds along the directions of an x axis, a y axis and a z axis under a machine body coordinate system are respectively set as u_w、v_w、w_wThree-dimensional wind speed u_w、v_w、w_wRespectively is

System noise $w\left(t\right)={[{\mathrm{\&delta;}}_{\stackrel{\&OverBar;}{L}},{\mathrm{\&delta;}}_M,{\mathrm{\&delta;}}_N,{\mathrm{\&delta;}}_D,{\mathrm{\&delta;}}_Y,{\mathrm{\&delta;}}_L,{\mathrm{\&delta;}}_{T_x},{\mathrm{\&delta;}}_{T_y},{\mathrm{\&delta;}}_{T_z},{\mathrm{\&delta;}}_{u_w},{\mathrm{\&delta;}}_{v_w},{\mathrm{\&delta;}}_{w_w},{\mathrm{\&delta;}}_{{\stackrel{.}{u}}_w},{\mathrm{\&delta;}}_{{\stackrel{.}{v}}_w},{\mathrm{\&delta;}}_{{\stackrel{.}{w}}_w},{\mathrm{\&delta;}}_m]}^{T^{\&prime;}}$ Respectively by external torque of aircraftM, N, aerodynamic force D, Y, L, engine thrust T_x、T_y、T_zWind speed u_w、v_w、w_wRate of change of wind speed

And total mass m, g is the acceleration of gravity, and the system noise figure array g (t) is defined as:

$G\left(t\right)=\left[\begin{array}{cccc}{0}_{2\&times;3}& 0_{2\&times;6}& 0_{2\&times;3}& 0_{2\&times;4}\\ A_{3\&times;3}& 0_{3\&times;6}& 0_{3\&times;3}& 0_{3\&times;4}\\ 0_{3\&times;3}& B_{3\&times;6}& C_{3\&times;3}& D_{3\&times;4}\end{array}\right]$

wherein, $A_{3\&times;3}=\left[\begin{array}{ccc}{c}_3& 0& c_4\\ 0& c_7& 0\\ c_4& 0& c_9\end{array}\right]$

$B_{3\&times;6}=\left[\begin{array}{cccccc}0& 0& \frac{-1}{\mathrm{mV}\cos \mathrm{\&beta;}}& \frac{-\sin \mathrm{\&alpha;}}{\mathrm{mV}\cos \mathrm{\&beta;}}& 0& \frac{\cos \mathrm{\&alpha;}}{\mathrm{mV}\cos \mathrm{\&beta;}}\\ 0& \frac{1}{\mathrm{mV}}& 0& \frac{-\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{\mathrm{mV}}& \frac{\cos \mathrm{\&beta;}}{\mathrm{mV}}& \frac{-\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{\mathrm{mV}}\\ \frac{-1}{m}& 0& 0& \frac{\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}}{m}& \frac{\sin \mathrm{\&beta;}}{m}& \frac{\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}}{m}\end{array}\right]$

$C_{3\&times;3}\left[\begin{array}{ccc}\frac{q\cos \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}& \frac{-p\cos \mathrm{\&alpha;}-r\sin \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}& \frac{q\sin \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}\\ \frac{-r\cos \mathrm{\&beta;}-q\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}& \frac{-r\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}+p\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}& \frac{p\cos \mathrm{\&beta;}+q\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}\\ -r\sin \mathrm{\&beta;}+q\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}& r\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}-p\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}& p\sin \mathrm{\&beta;}-q\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}\end{array}\right]$

$D_{3\&times;4}=\left[\begin{array}{cccc}\frac{\sin \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}& 0& \frac{-\cos \mathrm{\&alpha;}}{V\cos \mathrm{\&beta;}}& \frac{L+T_x\sin \mathrm{\&alpha;}-T_z\cos \mathrm{\&alpha;}}{m^{2}V\cos \mathrm{\&beta;}}\\ \frac{\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}& \frac{-\cos \mathrm{\&beta;}}{V}& \frac{\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{V}& \frac{-Y+T_x\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}-T_y\cos \mathrm{\&beta;}+T_z\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}}{m^{2}V}\\ -\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}& -\sin \mathrm{\&beta;}& -\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}& \frac{D-T_x\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}-T_y\sin \mathrm{\&beta;}-T_z\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}}{m^2}\end{array}\right]$

(b) establishment of extended Kalman filter measurement equation

Three-dimensional ground speed u obtained from GPS_g、v_g、w_gBy usingThe calculated V is used as the vacuum speed measurement and is combined with the pitch angle theta, the roll angle phi, the roll angle speed p, the pitch angle speed q, the yaw angle speed r and the three-dimensional linear acceleration a which are obtained by the strapdown inertial navigation system_x、a_y、a_zThe quantity of the total composition is measured Z, i.e. Z = [ theta, phi, p, q, r, a_x,a_y,a_z,V]^T' further establishing an extended Kalman filter measurement equation:

Z(t)=h[X(t),u(t),t]+v(t)

the measurement equation can be specifically expressed as:

$\left(\begin{array}{c}\mathrm{\&theta;}\\ \mathrm{\&phi;}\\ p\\ q\\ r\\ a_x\\ a_y\\ a_z\\ V\end{array}\right)=\left(\begin{array}{c}\mathrm{\&theta;}\\ \mathrm{\&phi;}\\ p\\ q\\ r\\ \frac{1}{m}(T_x+L\sin \mathrm{\&alpha;}-Y\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}-D\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;})\\ \frac{1}{m}(T_y+Y\cos \mathrm{\&beta;}-D\sin \mathrm{\&beta;})\\ \frac{1}{m}(T_z-L\cos \mathrm{\&alpha;}-Y\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-D\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;})\\ V\end{array}\right)+v\left(t\right)$

wherein the noise is measured

Caused by the measurement noise of the vacuum velocity calculated by the gyroscope and the linear accelerometer of the inertial device and the GPS;

(c) the state equation and the measurement equation are converted into a continuous linear interference equation:

converting the state equation in the step (a) and the measurement equation in the step (b) into a continuous linear interference equation form by using a Taylor series expansion method, namely f [ X (t), u (t), t)]And h [ X (t), u (t), t)]In-state optimal estimation value

And performing Taylor expansion, and omitting terms with more than two times to obtain a continuous linear interference equation:

$\left\{\begin{array}{c}\mathrm{\&delta;}\stackrel{\&CenterDot;}{X}\left(t\right)=F\left(t\right)\mathrm{\&delta;X}\left(t\right)+G\left(t\right)w\left(t\right)\\ \mathrm{\&delta;Z}\left(t\right)=H\left(t\right)\mathrm{\&delta;X}\left(t\right)+v\left(t\right)\end{array}\right.$

in the formula, $\mathrm{\&delta;X}\left(t\right)=[X\left(t\right)-\hat{X}\left(t\right)],$ $\mathrm{\&delta;Z}\left(t\right)=[Z\left(t\right)-\hat{Z}\left(t\right)],$ $F\left(t\right)={\frac{\&PartialD;f[X\left(t\right),u\left(t\right),t]}{\&PartialD;X\left(t\right)}|}_{X\left(t\right)=\hat{X}\left(t\right)},$ $H\left(t\right)={\frac{\&PartialD;h[X\left(t\right),u\left(t\right),t]}{\&PartialD;X\left(t\right)}|}_{X\left(t\right)=\hat{X}\left(t\right)},$ $\hat{Z}\left(t\right)=h[X\left(t\right),u\left(t\right),t]{|}_{X\left(t\right)=\hat{X}\left(t\right)};$

(d) discretizing a continuous linear disturbance equation:

sampling period Δ T = T_k+1-t_kDiscretizing the continuous linear disturbance equation obtained in the step (c) to obtain a discrete linear disturbance equation:

$\left\{\begin{array}{c}\mathrm{\&delta;}{X}_{k+1}={\mathrm{\&Phi;}}_{k+1,k}{\mathrm{\&delta;X}}_k+G_k{W}_k\\ \mathrm{\&delta;}{Z}_{k+1}=H_{k+1}\mathrm{\&delta;}{X}_{k+1}+V_{k+1}\end{array}\right.$

in the formula, ${\mathrm{\&delta;X}}_{k+1}=X_{k+1}-{\hat{X}}_{k+1/k},$ ${\mathrm{\&delta;Z}}_{k+1}=Z_{k+1}-{\hat{Z}}_{k+1},$ phi when sampling period DeltaT is small_k+1,k=I+F(t_k)·ΔT， $G_k=\mathrm{\&Delta;T}\&CenterDot;(I+\frac{F\left(t_k\right)}{2!}\&CenterDot;\mathrm{\&Delta;T})\&CenterDot;G\left(t_k\right),$ $F\left(t_k\right)={\frac{\&PartialD;f[X\left(t_k\right),u\left(t_k\right),t_k]}{\&PartialD;X\left(t_k\right)}|}_{X\left(t_k\right)={\hat{X}}_k},$

Is the last time instant t_kThe best estimate of the state of the moment in time,

for the current time instant, i.e. t_k+1One-step predicted value of the state quantity at the moment;

(e) extended kalman filter equation:

t is obtained according to the step (3)_kAerodynamic force D, Y, L at time, external torque

M, N, engine thrust T_x、T_y、T_zThe coefficient c of the moment equation set obtained in the step (7)₁、c₂、c₃、c₄、c₅、c₆、c₇、c₈、c₉T calculated in step (8)_kTotal mass m of the aircraft at the moment in time, and t_kThree-dimensional wind speed u of a moment_w、v_w、w_wAnd three dimensional rate of change of wind speed

Parameters in combination with t_kOptimal estimation of the state of a time

By using

Finding t_k+1One-step prediction value of time state quantity

t_k+1One-step prediction value of time state quantity

Variance matrix P of_k+1/kPassing through type

Solving is carried out; t is t_k+1Time of day filter gain matrix K_k+1Passing through type $K_{k+1}=P_{k+1/k}{H}_{k+1}^{T^{\&prime;}}{[H_{k+1}{P}_{k+1/k}{H}_{k+1}^{T^{\&prime;}}+R_{k+1}]}^{-1}$ Solving is carried out;

obtaining t according to the step (2)_k+1Three-dimensional ground speed u of time_g、v_g、w_gBy using

To obtain t_k+1Measuring the vacuum speed at the moment, and combining the t obtained in the step (1)_k+1Moment pitch angle theta, roll angle phi, roll angular velocity p, pitch angular velocity q, yaw angular velocity r, and three-dimensional linear acceleration a_x、a_y、a_zParameter composition t_k+1Total quantity of time measurement Z_k+1Combining the solved t_k+1One-step prediction value of time state quantityAnd a filter gain matrix K_k+1Using the formula $\mathrm{\&delta;}{\hat{X}}_{k+1}=K_{k+1}\{Z_{k+1}-h[{\hat{X}}_{k+1/k},u\left(t_{k+1}\right),t_{k+1}\left]\right\}$ Is carried out t_k+1Time state quantity improvement valueSolving;

will obtain t_k+1Time state quantity improvement value

And one-step prediction

Add, i.e. to

To obtain t_k+1Optimal estimation of the state of a time

Optimal estimated value

The error variance matrix of (2) can be represented by $P_{k+1}=(I-K_{k+1}{H}_{k+1})P_{k+1/k}{(I-K_{k+1}{H}_{k+1})}^{T^{\&prime;}}+K_{k+1}{R}_{k+1}{K}_{k+1}^{T^{\&prime;}}$ Solving is carried out, thereby realizing t under the flight state with large attack angle_k+1Accurately estimating the moment attack angle and the sideslip angle in real time;

in the extended Kalman filter, initial values of a pitch angle theta, a roll angle phi, a roll angular velocity p, a pitch angular velocity q and a yaw angular velocity r in state quantities are provided by a strapdown inertial navigation system inside a system, initial values of an attack angle alpha, a sideslip angle beta and a vacuum velocity V in the state quantities are provided by an atmospheric data system outside the system, and therefore an initial state is formed

Initial state estimation error variance matrix P₀Initial value Q of system noise variance array₀Measuring the initial value R of the noise variance matrix₀The initial values of the matrix are directly input from the outside of the system, and the noise matrix Q of the system_kInitial value Q of system noise variance matrix₀Determining and measuring a noise matrix R_k+1By measuring the initial value R of the noise variance matrix₀Determining, in addition, a three-dimensional wind speed u_w、v_w、w_wIs provided by an anemometer external to the system, three-dimensional wind speed rate of change

Are all taken to be zero.

4. The method according to claim 3, wherein the method for estimating the atmospheric angle of attack and the sideslip angle in the flight state with the large angle of attack comprises: the specific method for resolving the wind speed and the wind speed change rate in the step (10) is as follows:

t obtained according to step (1)_k+1The pitch angle theta, the roll angle phi and the yaw angle psi of the moment are obtained to obtain the body t_k+1Attitude transformation matrix S of time_θψφ：

$S_{\mathrm{\&theta;\&psi;\&phi;}}=\left[\begin{array}{ccc}\cos \mathrm{\&theta;}\cos \mathrm{\&psi;}& \cos \mathrm{\&theta;}\sin \mathrm{\&psi;}& -\sin \mathrm{\&theta;}\\ \sin \mathrm{\&theta;}\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}-\sin \mathrm{\&psi;}\cos \mathrm{\&phi;}& \sin \mathrm{\&theta;}\sin \mathrm{\&psi;}\sin \mathrm{\&phi;}+\cos \mathrm{\&psi;}\cos \mathrm{\&phi;}& \cos \mathrm{\&theta;}\sin \mathrm{\&phi;}\\ \sin \mathrm{\&theta;}\cos \mathrm{\&psi;}\cos \mathrm{\&phi;}+\sin \mathrm{\&psi;}\sin \mathrm{\&phi;}& \sin \mathrm{\&theta;}\sin \mathrm{\&psi;}\cos \mathrm{\&phi;}-\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}& \cos \mathrm{\&theta;}\cos \mathrm{\&phi;}\end{array}\right]\text{;}$

Combining t obtained in the step (2)_k+1Three-dimensional of time of dayGround speed u_g、v_g、w_gTo obtain t_k+1Three-dimensional ground speed u along x-axis, y-axis and z-axis directions under time machine body coordinate system_e,v_e,w_eComprises the following steps:

$\left[\begin{array}{c}{u}_e\\ v_e\\ w_e\end{array}\right]=S_{\mathrm{\&theta;\&psi;\&phi;}}\left[\begin{array}{c}{u}_g\\ v_g\\ w_g\end{array}\right];$

combining the t estimated in the step (9)_k+1The attack angle alpha, the sideslip angle beta and the vacuum speed V parameter of the moment are converted into a three-dimensional vacuum speed u along the directions of an x axis, a y axis and a z axis under a body coordinate system by using the following formula_b、v_b、w_bI.e. by

$\left[\begin{array}{c}{u}_b\\ v_b\\ w_b\end{array}\right]=\left[\begin{array}{c}V\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}\\ V\sin \mathrm{\&beta;}\\ V\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}\end{array}\right];$

Then t_k+1Three-dimensional wind speed u along x-axis, y-axis and z-axis directions under time machine body coordinate system_w、v_w、w_wThe calculation is made using the following formula:

$\left[\begin{array}{c}{u}_w\\ v_w\\ w_w\end{array}\right]=\left[\begin{array}{c}{u}_e\\ v_e\\ w_e\end{array}\right]-\left[\begin{array}{c}{u}_b\\ v_b\\ w_b\end{array}\right];$

using calculated t_k+1Three-dimensional wind speed u of a moment_w(t_k+1)、v_w(t_k+1)、w_w(t_k+1) Combined with t_kThree-dimensional wind speed u of a moment_w(t_k)、v_w(t_k)、w_w(t_k) By differentiating to find t_k+1Time three-dimensional wind speed rate of change

Namely:

$\left\{\begin{array}{c}{\stackrel{.}{u}}_w\left(t_{k+1}\right)=\frac{u_w\left(t_{k+1}\right)-u_w\left(t_k\right)}{\mathrm{\&Delta;T}}\\ {\stackrel{.}{v}}_w\left(t_{k+1}\right)=\frac{v_w\left(t_{k+1}\right)-v_w\left(t_k\right)}{\mathrm{\&Delta;T}}\\ {\stackrel{.}{w}}_w\left(t_{k+1}\right)=\frac{w_w\left(t_{k+1}\right)-w_w\left(t_k\right)}{\mathrm{\&Delta;T}}\end{array}\right..$

Images