High-precision off-line estimation method for aircraft airflow angle based on EM algorithm
1. An aircraft airflow angle high-precision off-line estimation method based on an EM algorithm is characterized by comprising the following steps:
step 1: carrying out stress analysis on the aircraft, and modeling the aircraft by utilizing a kinematic equation and a kinetic equation to obtain an adaptive state space model;
step 1-1: through a discrete body axis kinetic equation and a rotational kinematic equation, a state equation is obtained as follows:
rewriting formula (1) as:
xk+1＝f(xk，uk)+ωk (2)
in the formulaIs an unknown state quantity to be estimated, where θk，φkPitch and roll angles, u, respectivelyk，vk，wkRespectively is the projection of the airspeed in three directions of the body axis system; u. ofk＝[axk ayk azk pk qk rk]TIs a known input quantity, where pk，qk，rkRoll, pitch and yaw rates, axk，ayk，azkAcceleration values of the body axis system measured by the acceleration sensor in three directions are measured respectively; omegak＝[ωk(1) ωk(2) ωk(3) ωk(4) ωk(5)]TFor state noise, i.e. ω, following a Gaussian distribution of zero mean and unknown variancekN (0, Q), and Q is a positive definite matrix;
step 1-2: assuming that the aircraft is equipped with an airspeed sensor and a platform inertial navigation system, the measurement equation under windless conditions is established as equation (16):
the above formula is rewritten as:
zk＝h(xk，uk)+ζk (4)
whereinMeasuring the quantity at the moment k, wherein V is the airspeed of the aircraft; zetak＝[ζk(1) ζk(2) ζk(3)ζk(4) ζk(5) ζk(6)]TMeasurement noise, i.e.. zeta.kN (0, R), and R is a positive definite matrix;
the state space model of the nonlinear discrete system of the aircraft is formed by the formula (2) and the formula (4);
if the velocity component uk，vk，wkAs is known, the aircraft angle of attack α and sideslip angle β are calculated from equation (5):
step 2: adding different state noises and process noises by using the aircraft state space model established in the step 1 to obtain flight data in different noise environments, and generating simulation data for testing;
step 2-1: firstly, a group of state initial values are given, state quantities at other moments are obtained through a Runge Kutta method, and artificially set noise is added to obtain state quantities containing noise interference;
step 2-2: calculating a group of state quantities without noise interference by a Runge Kutta method, and calculating to obtain true values of an aircraft attack angle and a sideslip angle by the formula (5);
step 2-3: substituting the state quantity obtained in the step 2-2 into the formula (2), and adding the artificially set measured noise to obtain the observed quantity with noise;
step 2-4: changing the state noise and the measurement noise during each simulation, namely changing the Q and R to obtain simulation values under different noise levels;
and step 3: loading the simulation data generated in the step (2), operating an EM (effective noise) algorithm to obtain airflow angle estimation values under different noise backgrounds, comparing the estimation values with real values, and verifying the feasibility of the EM algorithm through Monte Carlo simulation;
step 3-1: applying the EM algorithm to a nonlinear discrete system state space model, and applying a state vector XN＝{x1，x2，…，xNAs unobservable data in EM, measure vector ZN＝{z1，z2，…，zNAs observable data in the EM,is unknown statistic;P1is a covariance matrix of the state quantities at the initial time,predicting a mean value for the state quantity at the initial moment;
by the EM algorithm, the cost function is as follows:
wherein the content of the first and second substances,is an estimate of the unknown statistical quantity,to give a measurement data ZNAnd current unknown statistic estimateLatent variable data XNConditional probability distribution of (1), pμ(XN，ZN) State vector X established for using mu as model parameterNAnd a measurement vector ZNA joint probability distribution of (a);
step 3-2: given a parameter μ, the log joint probability distribution of system states and observations is:
substituting formula (7) for formula (6) yields:
wherein
Step 3-3: the optimal estimation of unknown statistic is deduced by a maximized cost function, under the assumption of Gaussian distribution, a third-order spherical approximation algorithm is adopted to calculate the formulas (9) to (11), and the statistic mu has an analytic solution as follows:
P1＝P1∣N (13)
in the formula
i＝1，2，…4n
Wherein N is the total number of samples; xiiSigma points selected according to the sphere volume criterion; n is the dimension of the state vector; m is the dimension of the measurement;
and 4, step 4: and setting various parameters of the aircraft to obtain a corresponding estimated value of the airflow angle.
Background
The angle of attack and the angle of sideslip are two important flight state parameters of flight mechanics, and are two key parameters required by a flight control and navigation system. Due to the limited payload and budget, many aircraft systems, particularly drones, are not equipped with an airflow angle sensor, so that using data of other related sensors (such as GPS and INS) to estimate two airflow angles becomes an urgent problem to be solved.
In general, the airflow angle is directly measured by a vane sensor, a differential pressure sensor, a zero differential pressure sensor, and the like installed on an aircraft, but such sensors are expensive to use, and almost inevitably cause a large zero point deviation due to the influence of icing and local circulation currents related to the flight state. A tolerance deviation of the angle of attack will cause a large error in the calculation of the control law for the aircraft and may lead to an uncontrollable aircraft in certain angle of attack situations. In addition, the attack angle is also an important parameter required by an airplane flight control system and a stall warning, and the deviation of the sideslip angle can cause unnecessary energy loss of the airplane.
In the presence of system noise in an aircraft dynamic system, pneumatic parameter identification becomes a typical joint estimation problem of states and parameters. The state estimation and the parameter estimation are deeply coupled, mutually influenced and mutually causal; the state estimation error may cause a parameter identification error, which in turn may cause a state estimation error, thereby rendering the joint estimation problem more complex. This problem is further complicated by the need to extract a large amount of information from the experimental data, especially when the statistical properties of the system noise and the metrology noise are unknown. With the continuous development of the optimal estimation theory, the combination of an expectation maximization algorithm (EM algorithm) and a filtering/smoothing algorithm gradually becomes mature, the airflow angle estimation problem can be converted into a self-adaptive filtering problem, the numerical stability of the EM algorithm is good, and important research progress is made in the near term by applying the EM iterative optimization idea to realize the integrated processing of estimation and identification, so that the estimation of the process noise and the measurement noise of the system by adopting the EM algorithm is considered, the influence of the noise on the identification result can be effectively reduced, and the identification precision is greatly improved.
Disclosure of Invention
In order to overcome the defects of the prior art, the invention provides an aircraft airflow angle high-precision off-line estimation method based on an EM algorithm, which comprises the steps of firstly carrying out stress analysis on an aircraft, and modeling the aircraft by utilizing a kinematic equation and a kinetic equation to obtain an adaptive state space model; then, adding different state noises and process noises by using the established aircraft model to obtain flight data in different noise environments; then loading simulation data, operating an EM algorithm to obtain airflow angle estimation values under different noise backgrounds, comparing the estimation values with real values, and verifying the feasibility of the EM algorithm through Monte Carlo simulation; and finally, encapsulating the EM algorithm, so that a user can obtain a corresponding airflow angle estimation value only by setting various parameters of the aircraft. The method has higher estimation result precision and can meet the requirement of an aircraft tester on the precise airflow angle.
The technical scheme adopted by the invention for solving the technical problem comprises the following steps:
step 1: carrying out stress analysis on the aircraft, and modeling the aircraft by utilizing a kinematic equation and a kinetic equation to obtain an adaptive state space model;
step 1-1: through a discrete body axis kinetic equation and a rotational kinematic equation, a state equation is obtained as follows:
rewriting formula (1) as:
xk+1＝f(xk,uk)+ωk (2)
in the formulaIs an unknown state quantity to be estimated, where θk,φkPitch and roll angles, u, respectivelyk,vk,wkRespectively is the projection of the airspeed in three directions of the body axis system; u. ofk＝[axk ayk azk pk qk rk]TIs a known input quantity, where pk,qk,rkRoll, pitch and yaw rates, axk,ayk,azkAcceleration values of the body axis system measured by the acceleration sensor in three directions are measured respectively; omegak＝[ωk(1) ωk(2) ωk(3) ωk(4)ωk(5)]TFor state noise, i.e. ω, following a Gaussian distribution of zero mean and unknown variancekN (0, Q), and Q is a positive definite matrix;
step 1-2: assuming that the aircraft is equipped with an airspeed sensor and a platform inertial navigation system, the measurement equation under windless conditions is established as equation (16):
the above formula is rewritten as:
zk＝h(xk,uk)+ζk (4)
whereinMeasuring the quantity at the moment k, wherein V is the airspeed of the aircraft; zetak＝[ζk(1) ζk(2)ζk(3) ζk(4) ζk(5) ζk(6)]TMeasurement noise, i.e.. zeta.kN (0, R), and R is a positive definite matrix;
the state space model of the nonlinear discrete system of the aircraft is formed by the formula (2) and the formula (4);
if the velocity component uk,vk,wkAs is known, the aircraft angle of attack α and sideslip angle β are calculated from equation (5):
step 2: adding different state noises and process noises by using the aircraft state space model established in the step 1 to obtain flight data in different noise environments, and generating simulation data for testing;
step 2-1: firstly, a group of state initial values are given, state quantities at other moments are obtained through a Runge Kutta method, and artificially set noise is added to obtain state quantities containing noise interference;
step 2-2: calculating a group of state quantities without noise interference by a Runge Kutta method, and calculating to obtain true values of an aircraft attack angle and a sideslip angle by the formula (5);
step 2-3: substituting the state quantity obtained in the step 2-2 into the formula (2), and adding the artificially set measured noise to obtain the observed quantity with noise;
step 2-4: changing the state noise and the measurement noise during each simulation, namely changing the Q and R to obtain simulation values under different noise levels;
and step 3: loading the simulation data generated in the step (2), operating an EM (effective noise) algorithm to obtain airflow angle estimation values under different noise backgrounds, comparing the estimation values with real values, and verifying the feasibility of the EM algorithm through Monte Carlo simulation;
step 3-1: applying the EM algorithm to a nonlinear discrete system state space model, and applying a state vector XN＝{x1,x2,…,xNAs unobservable data in EM, measure vector ZN＝{z1,z2,…,zNAs observable data in the EM,is unknown statistic; p1Is a covariance matrix of the state quantities at the initial time,predicting a mean value for the state quantity at the initial moment;
by the EM algorithm, the cost function is as follows:
wherein the content of the first and second substances,is an estimate of the unknown statistical quantity,to give a measurement data ZNAnd current unknown statistic estimateLatent variable data XNConditional probability distribution of (1), pμ(XN,ZN) State vector X established for using mu as model parameterNAnd a measurement vector ZNA joint probability distribution of (a);
step 3-2: given a parameter μ, the log joint probability distribution of system states and observations is:
substituting formula (7) for formula (6) yields:
wherein
Step 3-3: the optimal estimation of unknown statistic is deduced by a maximized cost function, under the assumption of Gaussian distribution, a third-order spherical approximation algorithm is adopted to calculate the formulas (9) to (11), and the statistic mu has an analytic solution as follows:
P1＝P1∣N (13)
in the formula
i＝1,2,…4n
Wherein N is the total number of samples; xiiSigma points selected according to the sphere volume criterion; n is the dimension of the state vector; m is the dimension of the measurement;
and 4, step 4: and setting various parameters of the aircraft to obtain a corresponding estimated value of the airflow angle.
The invention has the following beneficial effects:
when the traditional airflow angle identification algorithm is used for estimating the attack angle and the sideslip angle, the variance statistical characteristics of state noise and measurement noise are needed, and inaccurate or even wrong noise assumption often greatly influences the accuracy of an estimation result, so that the final estimation result is influenced. According to the method, the statistical characteristic of the noise is estimated by using the expectation maximization algorithm, so that the reduction of the accuracy of the estimation result caused by the inaccuracy of the statistical characteristic of the noise is avoided. Compared with filtering, the smoothing method obtains higher-precision estimation by utilizing more measurement information, and in addition, the volume Kalman smoother has the advantages of high estimation precision, relatively small calculated amount and the like, so that the method is very suitable for the estimation problem of high-dimensional state quantity. As can be seen from the formulas (2) and (4), the algorithm does not use related information such as a pneumatic model and a pneumatic derivative when constructing the state space model of the nonlinear discrete system of the aircraft, so that the state space model is independent of the pneumatic layout of the aircraft, namely the estimation work of the airflow angle can be normally carried out under the special layout, the application range of the algorithm is expanded, and a foundation is laid for the identification of the parameters of the pneumatic model of the aircraft in the next step.
Drawings
FIG. 1 is an input volume change image used in generating simulation data according to an embodiment of the present invention.
Fig. 2 is a state quantity and measurement variation image obtained by flight simulation according to the embodiment of the present invention.
FIG. 3 shows the results of flow angle reconstruction using single simulation data according to an embodiment of the present invention.
Fig. 4 shows the results of 100 monte carlo simulations and the evaluation of each reconstruction result using the mean square error according to the embodiment of the present invention.
FIG. 5 shows the results of a comparison of an embodiment of the present invention with a conventional EKF algorithm.
FIG. 6 shows the noise identification result (in Q) according to the embodiment of the present inventionu,Qθ,RVAnd RθFor example).
FIG. 7 is a result of flow angle reconstruction using single real flight data for an embodiment of the present invention.
Detailed Description
The invention is further illustrated with reference to the following figures and examples.
The method utilizes the good numerical stability of the EM algorithm and the characteristic that a volume Kalman smoother (CKS) has high identification precision and good stability to the estimation problem in a high-dimensional state, so that a new airflow angle reconstruction algorithm based on the EM and CKS combined estimation (EM-CKS) is established and derived, the problems that a system equation in a dynamic system is complex in form, high in nonlinearity degree, always existing in noise and incapable of being eliminated and the like can be solved, and an effective tool is provided for estimating the attack angle and the sideslip angle of an aircraft.
The invention aims to provide a new idea for estimating an airflow angle for an aircraft, wherein a traditional estimation method is an Extended Kalman Filter (EKF) based on nonlinear kinematics and a measurement model, and the essential characteristic of the Kalman filter method is that a state space method is adopted to describe a system model, and the error between the system model and an actual dynamic system comprises system noise and measurement noise, namely the statistical characteristics of the process noise and the measurement noise are used as prior information to participate in estimation. These statistical characteristics are usually obtained by sensor ground tests or expert experience, and have a certain deviation, resulting in a low accuracy of the final estimation result. The method can well estimate the airflow angle under the condition of unknown noise, and the result precision is high.
An aircraft airflow angle high-precision off-line estimation method based on an EM algorithm comprises the following steps:
step 1: carrying out stress analysis on the aircraft, and modeling the aircraft by utilizing a kinematic equation and a kinetic equation to obtain an adaptive state space model;
step 1-1: through a discrete body axis kinetic equation and a rotational kinematic equation, a state equation is obtained as follows:
rewrite equation (1) to a concise form:
xk+1＝f(xk,uk)+ωk (2)
in the formulaIs an unknown state quantity to be estimated, where θk,φkPitch and roll angles, u, respectivelyk,vk,wkRespectively is the projection of the airspeed in three directions of the body axis system; u. ofk＝[axk ayk azk pk qk rk]TIs a known input quantity, where pk,qk,rkRoll, pitch and yaw rates, axk,ayk,azkAcceleration values of the body axis system measured by the acceleration sensor in three directions are measured respectively; omegak＝[ωk(1) ωk(2) ωk(3) ωk(4)ωk(5)]TFor state noise, i.e. ω, following a Gaussian distribution of zero mean and unknown variancekN (0, Q), and Q is a positive definite matrix;
step 1-2: assuming that the aircraft is equipped with an airspeed sensor and a platform inertial navigation system, the measurement equation under windless conditions is established as equation (16):
the above formula is rewritten as a concise form:
zk＝h(xk,uk)+ζk (4)
whereinMeasuring the quantity at the moment k, wherein V is the airspeed of the aircraft; zetak＝[ζk(1) ζk(2)ζk(3) ζk(4) ζk(5) ζk(6)]TMeasurement noise, i.e.. zeta.kN (0, R), and R is a positive definite matrix;
the state space model of the nonlinear discrete system of the aircraft is formed by the formula (2) and the formula (4);
if the velocity component uk,vk,wkAs is known, the aircraft angle of attack α and sideslip angle β are calculated from equation (5):
step 2: adding different state noises and process noises by using the aircraft state space model established in the step 1 to obtain flight data in different noise environments, and generating simulation data for testing;
step 2-1: firstly, a group of state initial values are given, state quantities at other moments are obtained through a Runge Kutta method, and artificially set noise is added to obtain state quantities containing noise interference;
step 2-2: calculating a group of state quantities without noise interference by a Runge Kutta method, and calculating to obtain truth values of an aircraft attack angle and a sideslip angle by a formula (5) as truth values for measuring the quality of an airflow angle reconstruction result obtained by the estimation method;
step 2-3: substituting the state quantity obtained in the step 2-2 into the formula (2), and adding the artificially set measured noise to obtain the observed quantity with noise;
step 2-4: changing the state noise and the measurement noise during each simulation, namely changing the Q and R to obtain simulation values under different noise levels;
as can be seen from the state equation and the measurement equation, this embodiment is performed in the absence of wind, with unknown system noise and measurement noise, and without an aircraft equipped with an airflow angle sensor;
and step 3: loading the simulation data generated in the step (2), operating an EM (effective noise) algorithm to obtain airflow angle estimation values under different noise backgrounds, comparing the estimation values with real values, and verifying the feasibility of the EM algorithm through Monte Carlo simulation;
step 3-1: applying the EM algorithm to a nonlinear discrete system state space model, and applying a state vector XN＝{x1,x2,…,xNAs unobservable data in EM, measure vector ZN＝{z1,z2,…,zNAs observable data in the EM,is unknown statistic;
by the EM algorithm, the cost function is as follows:
wherein the content of the first and second substances,is an estimated value of the unknown statistic;
step 3-2: given a parameter μ, the log joint probability distribution of system states and observations is:
substituting formula (7) for formula (6) yields:
wherein
Step 3-3: the optimal estimation of unknown statistic is deduced by a maximized cost function, under the assumption of Gaussian distribution, a third-order spherical approximation algorithm is adopted to calculate the formulas (9) to (11), and the statistic mu has an analytic solution as follows:
P1＝P1∣N (13)
in the formula
i＝1,2,…4n
Wherein N is the total number of samples; xiiSigma points selected according to the sphere volume criterion; n is the dimension of the state vector; m is the dimension of the measurement;
for the sake of brevity and clarity, the algorithm steps set forth herein are summarized as follows:
givingP1Initial values of Q and R;
executing the processes of volume Kalman filtering (CKF) and volume Kalman smoothing (CRTSS) to calculateAnd
executing EM algorithm and updatingP1Q and R;
fourthly, repeating the step one to the step three until a convergence condition is reached;
calculating values of an attack angle and a sideslip angle through the formula (5);
and 4, step 4: and setting various parameters of the aircraft to obtain a corresponding estimated value of the airflow angle.
Fig. 1 is a diagram showing an input quantity change image used for generating simulation data according to an embodiment of the present invention, and fig. 2 is a diagram showing a state quantity and quantity measurement change image obtained by performing flight simulation according to an embodiment of the present invention.
As can be seen from fig. 3, the single flight data generated randomly is used to verify the algorithm, so that ideal estimation results of the attack angle and the sideslip angle can be obtained. However, single trials are occasional and do not sufficiently illustrate the applicability of the CKS-EM algorithm to the problem of flow angle reconstruction. Therefore, in order to eliminate the contingency of a single experiment and increase the reliability and credibility of the algorithm, 100 Monte Carlo simulations are performed in the example of the algorithm, and the MSE (mean Square error) is used for evaluating the estimation result of the sensor error parameter at each time. As can be seen from fig. 4, in one hundred-time monte carlo simulation, the magnitude of the mean square error of each time is nine bits after the decimal point no matter the attack angle or the sideslip angle, which indicates that for any group of randomly generated airflow angles, the CKS-EM algorithm can reconstruct the attack angle and the sideslip angle very close to the true value, thereby not only verifying that an ideal airflow angle reconstruction value can be obtained by the algorithm, but also proving that the algorithm has certain universality. In addition, the method of airflow angle reconstruction using the CKS-EM algorithm described herein is somewhat generalizable, as it is true for any profile aircraft.
Fig. 5 shows that the system noise variance matrix and the measured noise variance matrix are respectively enlarged by 10 times and 20 times, substituted into the EKF algorithm for filtering to obtain a state quantity close to a true value, and then a reconstructed attack angle and a reconstructed sideslip angle are obtained. Repeating the above process for 100 times (i.e. performing 100 Monte Carlo simulations) to obtain the mean square error of each time, and comparing with the result obtained by the CKS-EM algorithm. FIG. 6 shows the noise identification result (in Q) according to the embodiment of the present inventionu,Qθ,RVAnd RθFor example).
It can be seen from fig. 7 that the present algorithm also has good estimation performance for real flight data.
