Aircraft airspeed angle of attack estimation method based on augmented generalized Kalman filter

Abstract

The present invention relates to the field of aircraft control technology, specifically to an aircraft airspeed angle of attack estimation method based on augmented generalized Kalman filtering, comprising the following steps: acquiring sensor data and correcting acceleration, calculating an initial estimate of the airspeed angle of attack, constructing an error matrix and scoring it, extracting a change trend adjustment parameter, calculating a state change coefficient, and fusing multi-source information to output a joint airspeed angle of attack quantity. In the present invention, an initial airspeed estimate is constructed by projecting and back-projecting the acceleration component onto the body coordinate system, combining the height difference to enhance physical matching, establishing a dynamic compensation mechanism based on slope scoring and trend change extraction, improving adaptability to non-stationary flight states, introducing pitch angular velocity, rudder deflection, and aerodynamic coefficient bias to make state evolution more structured, establishing a change coefficient through ratio calculation, strengthening the coupling between state prediction and maneuvering dynamics, weighted fusion of observation residuals and predicted values, and improving the convergence ability of the estimate to complex aerodynamic interference.

Description

Technical Field

[0001] The present invention relates to the technical field of aircraft control, and in particular to an aircraft airspeed angle of attack estimation method based on augmented generalized Kalman filtering. Background Art

[0002] The field of aircraft control technology encompasses improving aircraft stability and maneuverability, primarily involving the monitoring, estimation, and control of key parameters such as airspeed, attitude, heading, and angle of attack. The core of aircraft control involves precisely controlling the aircraft's flight state through various sensors, measurement devices, and control systems, ensuring safe flight in complex environments. This technical field is widely used in the design and operation of various aircraft, including civil aviation, military aircraft, and drones, encompassing multiple sub-areas such as autopilot, flight control systems, and flight state estimation. The key objectives of aircraft control technology are to enhance the autonomy, reliability, and safety of aircraft.

[0003] Among these, aircraft airspeed and angle of attack estimation methods utilize aircraft sensor information and control algorithms to estimate an aircraft's airspeed and angle of attack in real time. This technology specifically addresses the issue of accurate measurement of an aircraft's airspeed and angle of attack parameters during flight, particularly when direct measurement is difficult. By processing data from multiple aircraft sensors and incorporating mathematical models, estimation methods dynamically determine the aircraft's airspeed and angle of attack, providing real-time feedback to the aircraft's automatic control system. These methods generally achieve accurate airspeed and angle of attack estimation by fusing and processing sensor data and employing techniques such as Kalman filtering and least squares methods.

[0004] Relying on static sensor fusion and model estimation makes it difficult to characterize the dynamic nonlinear response during high-speed flight, resulting in a lag in airspeed and angle of attack estimates to sudden attitude changes. Time-domain data changes are not fully utilized, resulting in an inability of the estimation results to identify and adjust to changing trends. Inadequate structural state modeling leads to inaccurate predictions in multi-axis control or complex aerodynamic environments, affecting the consistency between estimated and actual flight states, causing deviations in control responses and reducing flight stability. Summary of the Invention

[0005] The purpose of the present invention is to solve the shortcomings of the prior art and propose an aircraft airspeed angle of attack estimation method based on augmented generalized Kalman filtering.

[0006] To achieve the above object, the present invention adopts the following technical solution: an aircraft airspeed angle of attack estimation method based on augmented generalized Kalman filtering, comprising the following steps:

[0007] S1: Obtain gyroscope angular velocity, three-axis acceleration, and altimeter data, project the three-axis acceleration components in the aircraft coordinate system, and then perform back-projection correction. Calculate an initial airspeed estimate based on the altitude difference, calculate an initial angle of attack estimate based on the three-axis acceleration components and the angle between them and gravity, and generate gyroscope state parameter values.

[0008] S2: Based on the gyroscope state parameter values, extract the airspeed estimate and angle of attack estimate within the estimation period, and the point-by-point difference between the airspeed estimate and angle of attack estimate within the current period, normalize them separately, and then combine them to construct an input matrix to generate a standardized error matrix value;

[0009] S3: Calling the standardized error matrix value, extracting time series row vectors in sequence, comparing the slope differences between adjacent time series row vectors, dividing the scoring level intervals, assigning scores to the time series row vectors, extracting the score change trend, selecting the slope increase value corresponding to the score change interval as a compensation factor, adjusting the airspeed and angle of attack in the gyroscope state parameter values, and generating an offset dynamic prediction result;

[0010] S4: Based on the offset dynamic prediction result, a state variable set is formed, the pitch angular velocity, the rudder deflection amount, and the lift and drag coefficient bias value are used as additional quantities, the slope of each state element is extracted and then calculated with the period interval ratio to generate a structural state change coefficient;

[0011] S5: Call the airspeed change value and angle of attack change value in the structural state change coefficient, combine the current observation projection residual and the predicted value, extract the relevant elements of the state error covariance matrix to set the weighted ratio, and fuse the data from each source in proportion to generate a joint airspeed angle of attack output.

[0012] As a further solution of the present invention, the standardized error matrix value includes a normalized airspeed error vector, a normalized angle of attack error vector, and an estimation period offset amplitude coefficient; the offset dynamic prediction result is specifically a corrected airspeed prediction value, a corrected angle of attack prediction value, and an airspeed and angle of attack change trend fitting factor; the structural state change coefficient includes an airspeed change slope coefficient, an angle of attack change slope coefficient, and a state response time ratio; the joint airspeed and angle of attack output includes an airspeed estimation fusion result, an angle of attack estimation fusion result, and a residual weighted fusion coefficient.

[0013] As a further solution of the present invention, the step of obtaining the standardized error matrix value is specifically as follows:

[0014] S201: Based on the gyroscope state parameter value, obtain the airspeed estimate and the angle of attack estimate within the estimation period, calculate the point-by-point difference between the current period airspeed value and the estimated period airspeed estimate, and combine the point-by-point difference between the current period angle of attack value and the estimated period angle of attack estimate to generate a period difference sequence matrix;

[0015] S202: calling the periodic difference sequence matrix, dividing the multi-periodic point difference values by the corresponding maximum difference value for normalization, combining the normalized airspeed difference values and the angle of attack difference values to construct a matrix, and obtaining a normalized combined input matrix;

[0016] S203: Based on the normalized combined input matrix, the normalized airspeed difference and angle of attack difference of multiple period points are called, and the difference with the root mean square value and the average value of the corresponding channel is calculated using the formula:

[0017] ;

[0018] The standardized error values of multiple period points are obtained by operation and combined to generate the standardized error matrix value;

[0019] in, , represents the normalized airspeed difference at the qth period point, represents the normalized angle of attack difference at the qth period point, represents the RMS value of the angle of attack difference sequence of the p-th channel in the current estimation cycle, Represents the average value of the normalized combined input matrix of the p-th channel in the current estimation cycle.

[0020] As a further solution of the present invention, the step of obtaining the offset dynamic prediction result is specifically as follows:

[0021] S301: calling the standardized error matrix value, extracting adjacent time series row vectors according to the row vector index order, calculating the slope difference of corresponding elements between adjacent time series row vectors, taking the absolute value and then the average, and generating a slope difference value;

[0022] S302: Based on the slope difference values, divide the scoring level intervals according to a preset threshold range, compare each slope difference value with the boundary values of multiple intervals step by step, store the scoring values corresponding to the matching intervals into an array, and generate a scoring vector;

[0023] S303: Call the rating vector to calculate the incremental difference between adjacent rating values using the formula:

[0024] ;

[0025] The compensation factor is obtained by calculation and linearly superimposed with the initial airspeed nominal value and the angle of attack baseline parameter to generate the offset dynamic prediction result;

[0026] in, represents the compensation factor of the qth time series row vector at time s, Represents the slope difference value of adjacent time series row vectors, represents the rating vector, Represents the difference between the mth adjacent rating changes, represents the length of the rating vector, represents the time interval between adjacent time series row vectors, represents the angle of attack adjustment in the previous iteration cycle, Represents the upper threshold of the angle of attack

[0027] As a further solution of the present invention, the step of obtaining the structural state variation coefficient is specifically as follows:

[0028] S401: Based on the offset dynamic prediction result, the pitch angular velocity, the control surface deflection, and the lift and drag coefficient offset are called, and the three are aligned with the lateral displacement data and longitudinal acceleration data of the gyroscope state parameters in time series. The three are combined into a multidimensional set of dynamic parameters and basic state parameters to generate an extended set of state variables;

[0029] S402: extracting the time-domain variation slope of the pitch angular velocity, the time-domain variation slope of the rudder deflection, and the time-domain variation slope of the lift and drag coefficient bias value from the state variable extended set, performing linear fitting on the time-domain variation slopes of the three parameters using the least squares method, and calculating the absolute value of the variation rate to construct a state element slope set;

[0030] S403: Calling the time-domain variation slope of the pitch angular velocity, the time-domain variation slope of the rudder deflection, and the time-domain variation slope of the lift and drag coefficient bias value in the state element slope set, and calculating the ratio of the absolute value of the multiple slopes to the period interval parameter using the formula:

[0031] ;

[0032] By calculating the dynamic trend principal component and the rudder coupling compensation term, the structural state variation coefficient parameters are generated by superposition;

[0033] in, represents the structural state variation coefficient parameter, Represents the slope value of the i-th element in the state element slope set, represents the period interval parameter, Represents the minimum parameter of the lift-drag coefficient bias value period, Represents the maximum value parameter of the rudder deflection, represents the dynamic coupling factor parameter of the rudder surface, represents the periodic average value parameter of the lift-drag coefficient bias value, Represents the lift-drag stability compensation parameter

[0034] As a further solution of the present invention, the step of obtaining the combined airspeed angle of attack output is specifically as follows:

[0035] S501: calling the airspeed change value in the structural state change coefficient, extracting the transverse covariance element in the state error covariance matrix, calculating the ratio of the observation projection residual to the predicted value, setting the weighted proportional coefficient of the transverse covariance element and the longitudinal covariance element, and generating the covariance matrix weighted proportional coefficient;

[0036] S502: The weighted proportional coefficient of the covariance matrix is called, and the sum of the absolute differences between the observed projection residual and the predicted value is combined to calculate the square root of the square of the airspeed change increment and the angle of attack change increment using the formula:

[0037] ;

[0038] Calculate the weighted fusion factor, integrate the vertical elements of the covariance matrix, and generate the multi-source data fusion weight;

[0039] in, represents the m-th dimension weighted fusion factor, represents the m-th dimension covariance weighted proportional coefficient, represents the qth group of observation projection residuals, represents the predicted value of group q, represents the time series airspeed change increment, represents the increment of the angle of attack change in the spatial sequence, represents the horizontal elements in the covariance matrix, represents the vertical element in the covariance matrix, is the denominator protection constant;

[0040] S503: Calling the airspeed change value and the angle of attack change value, calculating the weighted summation result according to the multi-source data fusion weight, performing error threshold boundary truncation processing on the fusion result, and generating a joint airspeed angle of attack output.

[0041] As a further embodiment of the present invention, the method further comprises:

[0042] S1: Obtain gyroscope angular velocity, three-axis acceleration, and altimeter data, project the acceleration components in the aircraft coordinate system, and then perform back-projection correction. Calculate the initial airspeed estimate based on the altitude difference, calculate the initial angle of attack estimate based on the acceleration components and the angle between them and gravity, and generate gyroscope state parameter values.

[0043] The gyroscope state parameter values specifically include the body attitude angular rate, acceleration correction residual, and initial angle of attack estimation value.

[0044] As a further solution of the present invention, the step of obtaining the gyroscope state parameter value is specifically as follows:

[0045] S101: Based on the gyroscope angular velocity, three-axis acceleration, and altimeter data, the linear acceleration is obtained and projected into the aircraft coordinate system. Combined with the attitude angle, it is back-projected into the inertial system. The velocity change is calculated based on the altitude difference to generate an initial airspeed estimate.

[0046] S102: Call the initial airspeed estimate and the longitudinal acceleration component, calculate the angle with the gravity direction, and determine the incident direction and attitude offset using the formula:

[0047] ;

[0048] The calculation obtains the estimated value of the longitudinal incident angle, obtains the deviation trend based on the difference between it and the attitude angle, obtains the initial estimate of the angle of attack, and generates the gyroscope state parameter value;

[0049] in, represents the initial estimate of the angle of attack, Represents the acceleration component in the Z-axis direction in the body coordinate system, represents the acceleration due to gravity, represents the pitch angle, represents the initial estimate of airspeed, Respectively represent the acceleration components of the X-axis and Y-axis directions in the body coordinate system, Represents the altitude values obtained by the altimeter at the current moment and the previous moment respectively.

[0050] Compared with the prior art, the advantages and positive effects of the present invention are:

[0051] In the present invention, the acceleration components are projected and back-projected into the body coordinate system, and the initial airspeed estimate is constructed in combination with the altitude difference, thereby enhancing the physical matching of the original estimate. The normalization of the periodic difference makes the data input more consistent, and the slope score and trend change extraction establish a dynamic compensation mechanism to improve the adaptability to non-stationary flight conditions. The introduction of pitch angular velocity, rudder deflection and aerodynamic coefficient bias makes the state evolution more structured, and the variation coefficient is established through ratio calculation to strengthen the coupling between state prediction and maneuvering dynamics. The observation residual and the predicted value are weighted and fused to improve the convergence ability of the estimation to complex aerodynamic interference. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 It is a schematic diagram of the main steps of the present invention;

[0053] Figure 2 Flowchart of the steps for obtaining the state parameter value of the gyroscope according to the present invention;

[0054] Figure 3 Flowchart of the steps for obtaining the standardized error matrix value of the present invention;

[0055] Figure 4Flowchart of the steps for obtaining the offset dynamic prediction results of the present invention;

[0056] Figure 5 Flowchart of the steps for obtaining the structural state variation coefficient of the present invention;

[0057] Figure 6 The figure is a flow chart of the steps for obtaining the combined airspeed angle of attack output of the present invention. DETAILED DESCRIPTION

[0058] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0059] In the description of the present invention, it should be understood that the terms "length," "width," "up," "down," "front," "back," "left," "right," "vertical," "horizontal," "top," "bottom," "inside," "outside," and the like, indicating positions or relationships, are based on the positions or relationships shown in the accompanying drawings and are intended only to facilitate the description of the present invention and simplify the description. They do not indicate or imply that the devices or elements referred to must have a specific orientation, be constructed, or operate in a specific orientation. Therefore, they should not be construed as limiting the present invention. Furthermore, in the description of the present invention, "plurality" means two or more, unless otherwise expressly and specifically defined.

[0060] Example 1

[0061] See also Figure 1 The present invention provides a technical solution: an aircraft airspeed angle of attack estimation method based on augmented generalized Kalman filtering, comprising the following steps:

[0062] S1: Obtain gyroscope angular velocity, three-axis acceleration, and altimeter data, project the three-axis acceleration components in the aircraft coordinate system, and then perform back-projection correction. Calculate an initial airspeed estimate based on the altitude difference, calculate an initial angle of attack estimate based on the three-axis acceleration components and the angle between them and gravity, and generate gyroscope state parameter values.

[0063] S2: Based on the gyroscope state parameter values, extract the airspeed estimate and angle of attack estimate within the estimation period, and the point-by-point difference between the airspeed estimate and angle of attack estimate within the current period, normalize them separately, and then combine them to construct an input matrix to generate a standardized error matrix value;

[0064] S3: Calling the standardized error matrix value, extracting time series row vectors in sequence, comparing the slope differences between adjacent time series row vectors, dividing the scoring level intervals, assigning scores to the time series row vectors, extracting the score change trend, selecting the slope increase value corresponding to the score change interval as a compensation factor, adjusting the airspeed and angle of attack in the gyroscope state parameter values, and generating an offset dynamic prediction result;

[0065] S4: Based on the offset dynamic prediction result, a state variable set is formed, the pitch angular velocity, the rudder deflection amount, and the lift and drag coefficient bias value are used as additional quantities, the slope of each state element is extracted and then calculated with the period interval ratio to generate a structural state change coefficient;

[0066] S5: Call the airspeed change value and angle of attack change value in the structural state change coefficient, combine the current observation projection residual and the predicted value, extract the relevant elements of the state error covariance matrix to set the weighted ratio, and fuse the data from each source in proportion to generate a joint airspeed angle of attack output.

[0067] The gyroscope state parameter values are specifically the body attitude angular rate, acceleration correction residual, and initial angle of attack estimate. The standardized error matrix values include the normalized airspeed error vector, the normalized angle of attack error vector, and the estimation period offset amplitude coefficient. The offset dynamic prediction results are specifically the corrected airspeed prediction value, the corrected angle of attack prediction value, and the airspeed and angle of attack change trend fitting factor. The structural state change coefficients include the airspeed change slope coefficient, the angle of attack change slope coefficient, and the state response time ratio. The joint airspeed and angle of attack output includes the airspeed estimation fusion result, the angle of attack estimation fusion result, and the residual weighted fusion coefficient.

[0068] See also Figure 2 , the steps for obtaining the gyroscope state parameter value are as follows:

[0069] S101: Based on the gyroscope angular velocity, three-axis acceleration, and altimeter data, the linear acceleration is obtained and projected into the aircraft coordinate system. Combined with the attitude angle, it is back-projected into the inertial system. The velocity change is calculated based on the altitude difference to generate an initial airspeed estimate.

[0070] Based on gyroscope angular velocity 、 、 (measured by MEMS gyroscope, sampling frequency 100Hz), three-axis accelerometer raw data 、 、 (After low-pass filtering), altimeter data , do the following:

[0071] 1. Remove gravity: by attitude angle 、 (obtained by gyroscope integration and Kalman filter fusion), calculate the components of gravity on each axis 、 、 , and get the linear acceleration 、 、 ;

[0072] 2. Projection to inertial system: through rotation matrix (calculated by attitude angle), Projected as inertial acceleration ,in ;

[0073] 3. Integral speed change: combined with height difference (Time interval 1s), calculate vertical velocity , initial estimate of airspeed (corrected to 25m / s after smoothing filtering).

[0074] S102: Call the initial airspeed estimate and the longitudinal acceleration component, and calculate the angle with the gravity direction to determine the incident direction and attitude offset using the formula:

[0075] ;

[0076] The calculation obtains the estimated value of the longitudinal incident angle, obtains the deviation trend based on the difference between it and the attitude angle, obtains the initial estimate of the angle of attack, and generates the gyroscope state parameter value;

[0077] in, represents the initial estimate of the angle of attack, Represents the acceleration component in the Z-axis direction in the body coordinate system, represents the acceleration due to gravity, represents the pitch angle, represents the initial estimate of airspeed, Respectively represent the acceleration components of the X-axis and Y-axis directions in the body coordinate system, Represent the altitude values obtained by the altimeter at the current moment and the previous moment respectively;

[0078] During the calculation of the initial angle of attack, the initial airspeed estimate is first called , collect the Z-axis acceleration component in the body coordinate system (Measured by a three-axis accelerometer, the original data is , removing the gravity component Then get the pitch angle (Calculated by integrating the angular velocity of the gyroscope, the angular velocity , integration time , obtained after error compensation), gravitational acceleration , the X-axis and Y-axis acceleration components are 、 (obtained after coordinate projection), current altitude , the height at the previous moment (Measured by barometric altimeter, sampling interval is ).

[0079] The first term in the formula The calculation process is:

[0080] 1. Calculation ;

[0081] 2. Calculation ;

[0082] 3. Divide by get .

[0083] The second term in the formula The calculation process is:

[0084] 1. Calculation ;

[0085] 2. Calculation .

[0086] Add the two items and multiply them :

[0087] 1. Calculation ;

[0088] 2. Multiply get ;

[0089] 3. Calculate the square root of the height difference ;

[0090] 4. Final Result (about 31.5°).

[0091] Table 1: Parameters related to initial airspeed estimation

[0092]

[0093] As shown in Table 1, the initial airspeed estimate and associated parameters are obtained through sensor data and coordinate transformation. The formula quantifies the initial angle of attack estimate by combining the proportional relationship between longitudinal acceleration and lateral acceleration with the altitude change rate. This result indicates that the estimated angle of attack exceeds the range of normal flight conditions (typical value <15°) and needs to be corrected by combining the attitude angle difference.

[0094] See also Figure 3 , the specific steps for obtaining the standardized error matrix value are:

[0095] S201: Based on the gyroscope state parameter values, obtain the airspeed estimate and the angle of attack estimate within the estimation period, calculate the point-by-point difference between the current period airspeed value and the estimated period airspeed estimate, and combine the point-by-point difference between the current period angle of attack value and the estimated period angle of attack estimate to generate a period difference sequence matrix;

[0096] Based on gyroscope state parameter values , set the estimation period (Contains 5 data points, sampling interval 1s), call the current cycle airspeed value (measured by pitot tube and smoothed by Kalman filter), estimated periodic airspeed estimate (Through weighted prediction of state parameters and historical data), the angle of attack value of the current cycle , estimated period angle of attack , do the following:

[0097] 1. Calculate the point-by-point difference in airspeed ( to 5), we get ;

[0098] 2. Calculate the point-by-point difference in angle of attack ,get ;

[0099] 3. Arrange the difference sequence into a matrix form according to the periodic points:

[0100] ;

[0101] Table 2: Gyroscope status parameters and related data

[0102]

[0103] As shown in Table 2, the airspeed and angle of attack differences are generated by comparing the measured values with the predicted values point by point, and the matrix dimension is consistent with the number of periodic points.

[0104] S202: calling the periodic difference sequence matrix, dividing the multi-periodic point difference values by the corresponding maximum difference value for normalization, combining the normalized airspeed difference values and the angle of attack difference values to construct a matrix, and obtaining a normalized combined input matrix;

[0105] Call the periodic difference sequence matrix to extract the airspeed difference column , angle of attack difference column , do the following:

[0106] 1. Calculate the maximum airspeed difference , the maximum value of the angle of attack difference ;

[0107] 2. Normalized airspeed difference: ,get ;

[0108] 3. Normalized angle of attack difference: ,get ;

[0109] 4. Combine into a normalized input matrix:

[0110] ;

[0111] S203: Based on the normalized combined input matrix, call the normalized airspeed difference and angle of attack difference of multiple period points, and calculate the difference with the corresponding channel root mean square value and average value using the formula:

[0112] ;

[0113] The standardized error values of multiple period points are obtained by operation and combined to generate the standardized error matrix value;

[0114] in, represents the qth cycle.

[0115] Call the normalized combined input matrix and set the channel (airspeed difference) and (angle of attack difference), calculate the root mean square value of each channel and the average :

[0116] 1. For channel 1 (airspeed difference):

[0117] RMS value ;

[0118] average value ;

[0119] 2. For channel 2 (angle of attack difference):

[0120] RMS value ;

[0121] average value ;

[0122] Substitute the formula to calculate the standardized error value (by For example):

[0123] 1. , (normalized difference of channel 1);

[0124] 2. ;

[0125] 3. ,Right now ;

[0126] Similarly, calculate the error values of all periodic points and generate the standardized error matrix:

[0127] ;

[0128] The formula is beneficial because it dynamically adjusts the normalized difference weights by introducing channel root mean square values and average values, eliminating the interference of different dimensional parameters on the error assessment. The results show that the normalized error matrix can quantify the fluctuation characteristics of airspeed and angle of attack estimates within a cycle, providing an input basis for subsequent control compensation.

[0129] See also Figure 4 , the specific steps for obtaining the offset dynamic prediction results are:

[0130] S301: Call the standardized error matrix value, extract adjacent time series row vectors according to the row vector index order, calculate the slope difference of corresponding elements between adjacent time series row vectors, take the absolute value and then take the average to generate the slope difference value;

[0131] Call the standardized error matrix value and set the adjacent time series row vector index To 4 (the number of matrix rows is 5, and the interval between adjacent rows is 1 second), the error matrix is extracted as follows:

[0132] ;

[0133] Do the following:

[0134] 1. Calculate the slope difference of corresponding elements between adjacent row vectors:

[0135] right and , the slope difference of channel 1 (airspeed error) is , the slope difference of the second channel (angle of attack error) is , after taking the absolute value, the mean is ;

[0136] right and , the slope difference of channel 1 is , the second channel is , the mean is ;

[0137] right and , the slope difference of channel 1 is , the second channel is , the mean is ;

[0138] right and , the slope difference of channel 1 is , the second channel is , the mean is ;

[0139] 2. Generate a slope difference value sequence: .

[0140] Table 3: Example of a standardized error matrix

[0141]

[0142] As shown in Table 3, the slope differences of adjacent row vectors are calculated by channel-by-channel difference, and the mean reflects the intensity of time series fluctuations.

[0143] S302: Based on the slope difference values, the scoring level intervals are divided according to the preset threshold range, each slope difference value is compared with the boundary values of multiple intervals step by step, and the scoring values corresponding to the matching intervals are stored in an array to generate a scoring vector;

[0144] Based on the slope difference value series , preset rating range:

[0145] Low volatility (0≤difference<0.1): score 1;

[0146] Medium fluctuation (0.1≤difference<0.2): score 2;

[0147] High volatility (difference value ≥ 0.2): score 3;

[0148] Perform a level-by-level comparison:

[0149] 1. A difference of 0.1 is in the medium fluctuation range and is rated 2.

[0150] 2. A difference of 0.15 is in the medium fluctuation range and is rated 2;

[0151] 3. A difference of 0.3 is in the high volatility range and is rated 3;

[0152] 4. A difference of 0.25 is considered high volatility and is rated 3.

[0153] Generate a rating vector: .

[0154] Table 4: Matching table of scoring level intervals and difference values

[0155]

[0156] As shown in Table 4, the scoring level is achieved by comparing the interval boundary values step by step, and the difference values are mapped one-to-one with the scoring results.

[0157] S303: Call the rating vector and calculate the incremental difference between adjacent rating values using the formula:

[0158] ;

[0159] The compensation factor is obtained by calculation and linearly superimposed with the initial airspeed nominal value and the angle of attack baseline parameter to generate the offset dynamic prediction result;

[0160] in, represents the compensation factor of the qth time series row vector at time s, Represents the slope difference value of adjacent time series row vectors, represents the rating vector, Represents the difference between the mth adjacent rating changes, represents the length of the rating vector, represents the time interval between adjacent time series row vectors, represents the angle of attack adjustment in the previous iteration cycle, Represents the upper threshold of angle of attack.

[0161] Calling the score vector , the angle of attack adjustment amount of the previous iteration cycle , upper threshold of angle of attack (Set according to the aerodynamic performance of the aircraft), time interval , calculate the incremental difference of adjacent scores :

[0162] 1. ;

[0163] 2. ;

[0164] 3. ;

[0165] Substitute into the formula to calculate the compensation factor (in For example):

[0166] ;

[0167] Similarly, calculate the compensation factors for other periodic points and generate the results: , the compensation factor is added to the initial airspeed nominal value , angle of attack baseline Linear superposition (weight 0.1) gives the predicted offset:

[0168] Airspeed offset: ;

[0169] Angle of attack offset: .

[0170] The formula's benefit lies in its ability to dynamically adjust the compensation factor to account for time sensitivity and aerodynamic constraints by integrating the score increment with historical angle of attack adjustments. The results show that the compensation factor adjusts proportionally with the score level and angle of attack threshold, effectively mitigating overshoot risk.

[0171] See also Figure 5 , the steps for obtaining the structural state variation coefficient are as follows:

[0172] S401: Based on the offset dynamic prediction results, the pitch angular velocity, rudder deflection, and lift-drag coefficient offset values are called and aligned with the lateral displacement data and longitudinal acceleration data of the gyroscope state parameters in time series. These are then merged into a multidimensional set of dynamic parameters and basic state parameters to generate an extended set of state variables.

[0173] Based on the dynamic prediction results of the deviation (airspeed deviation , angle of attack deviation ), call the pitch angular velocity parameters (measured by gyroscope), rudder deflection parameter (Feedback from the servo), lift and drag coefficient offset value (calibrated by wind tunnel test), lateral displacement data (GPS measurement), longitudinal acceleration data (Accelerometer acquisition), aligned by time series (sampling interval ), merged into a multidimensional set:

[0174] ;

[0175] Table 5: Example table of multi-dimensional state variable extension set

[0176]

[0177] As shown in Table 5, the multidimensional set integrates dynamic prediction parameters and basic state parameters through timestamp alignment.

[0178] S402: Extracting the time-domain variation slope of the pitch angular velocity, the time-domain variation slope of the rudder deflection, and the time-domain variation slope of the lift and drag coefficient bias from the extended set of state variables, performing linear fitting on the time-domain variation slopes of the three parameters using the least squares method, and calculating the absolute value of the variation rate to construct a state element slope set;

[0179] Extract the pitch rate, rudder deflection, and lift and drag offset parameters from the extended set, and perform a least squares linear fit on the time series data (4 period points) of each parameter:

[0180] 1. Pitch rate slope:

[0181] Data points: , ;

[0182] Slope calculation: ;

[0183] 2. Rudder deflection slope:

[0184] Data points: ;

[0185] Slope calculation: ;

[0186] 3. Lift-drag bias slope:

[0187] Data points: ;

[0188] Slope calculation: ;

[0189] Construct a set of state element slopes: .

[0190] S403: Call the time-domain variation slope of the pitch angular velocity, the time-domain variation slope of the rudder deflection, and the time-domain variation slope of the lift and drag coefficient bias value in the state element slope set, and calculate the ratio of the absolute value of the multiple slopes to the period interval parameter using the formula:

[0191] ;

[0192] By calculating the dynamic trend principal component and the rudder coupling compensation term, the structural state variation coefficient parameters are generated by superposition;

[0193] in, represents the structural state variation coefficient parameter, Represents the slope value of the i-th element in the state element slope set, represents the period interval parameter, Represents the minimum parameter of the lift-drag coefficient bias value period, Represents the maximum value parameter of the rudder deflection, represents the dynamic coupling factor parameter of the rudder surface, represents the periodic average value parameter of the lift-drag coefficient bias value, Represents the lift-drag stability compensation parameter.

[0194] Call Slope Collection , period interval parameter , the minimum value of lift-drag coefficient during the bias period (Extracted from the fourth row of Table 1), the maximum deflection of the rudder surface , dynamic coupling factor of the rudder surface (calibrated according to the rudder efficiency curve), the periodic average value of the lift and drag coefficient , stability compensation (To prevent the denominator from being zero), substitute into the formula:

[0195] ;

[0196] Step-by-step calculation:

[0197] 1. Calculate the sum of squared slopes: ;

[0198] 2. Take the square root and multiply : ;

[0199] 3. Denominator: ;

[0200] 4. The first result: ;

[0201] 5. The second molecule: ;

[0202] 6. The second denominator: ;

[0203] 7. Second result: ;

[0204] 8. Final .

[0205] The formula is beneficial in that it quantifies the combined effects of aerodynamics and control by integrating the principal component of the dynamic trend with the rudder coupling effect. The results indicate that the coefficient of structural state variation is significantly dominated by the rudder deflection, necessitating prioritization of rudder compensation.

[0206] See also Figure 6 The specific steps for obtaining the combined airspeed angle of attack output are as follows:

[0207] S501: Call the airspeed change value in the structural state change coefficient, extract the transverse covariance element in the state error covariance matrix, calculate the ratio of the observation projection residual to the predicted value, set the weighted proportional coefficient of the transverse covariance element and the longitudinal covariance element, and generate the covariance matrix weighted proportional coefficient;

[0208] Call structure state change coefficient Airspeed change in , extract the horizontal covariance elements in the state error covariance matrix (corresponding to lateral displacement error), longitudinal covariance element (corresponding to longitudinal acceleration error), observation projection residual (GPS measured and predicted difference), predicted value , do the following:

[0209] 1. Calculate the ratio of observed residuals to predicted values: ;

[0210] 2. Set the horizontal and vertical covariance weighted proportional coefficients ;

[0211] 3. Generate a set of covariance matrix weighted proportional coefficients: .

[0212] Table 6: Example table of state error covariance matrix

[0213]

[0214] As shown in Table 6, the weighted proportional coefficient is directly calculated by the ratio of covariance elements.

[0215] S502: Call the covariance matrix weighted proportional coefficient, combine the absolute difference between the observed projection residual and the predicted value, and calculate the square root of the square of the airspeed change increment and the angle of attack change increment using the formula:

[0216] ;

[0217] Calculate the weighted fusion factor, integrate the vertical elements of the covariance matrix, and generate the multi-source data fusion weight;

[0218] in, represents the m-th dimension weighted fusion factor, represents the m-th dimension covariance weighted proportional coefficient, represents the qth group of observation projection residuals, represents the predicted value of group q, represents the time series airspeed change increment, represents the increment of the angle of attack change in the spatial sequence, represents the horizontal elements in the covariance matrix, represents the vertical element in the covariance matrix, is the denominator protection constant;

[0219] Call weighted proportional coefficient , the sum of the absolute differences between the observed residuals and the predicted values , airspeed change increment , angle of attack change increment (Convert to radians: ), denominator protection constant , substitute into the formula to calculate the weighted fusion factor (in For example):

[0220] ;

[0221] The same calculation , , generate multi-source data fusion weights: .

[0222] S503: Call the airspeed change value and the angle of attack change value, calculate the weighted summation result according to the multi-source data fusion weight, perform error threshold boundary truncation processing on the fusion result, and generate a joint airspeed angle of attack output.

[0223] Calling airspeed change value , angle of attack change value , fusion weight , perform a weighted sum:

[0224] ;

[0225] Set the error threshold boundary to (Based on the nominal airspeed of 25m / s, the threshold range is 23.75-26.25m / s), truncation processing:

[0226] If the combined output If it is within the threshold, it is output directly;

[0227] If exceeded, it will be truncated to the boundary value.

[0228] The formula is beneficial in that it balances the reliability of horizontal and vertical data through the covariance ratio and dynamic increment fusion. The results show that the weighting factor is dominated by the covariance ratio, ensuring that the weight of high-error dimensions is attenuated.

[0229] The above are merely preferred embodiments of the present invention and do not limit the present invention in any other form. Any technician familiar with the profession may use the technical content disclosed above to change or modify it into an equivalent embodiment with equivalent changes and apply it to other fields. However, any simple modification, equivalent change and modification made to the above embodiment based on the technical essence of the present invention without departing from the content of the technical solution of the present invention shall still fall within the scope of protection of the technical solution of the present invention.

Claims

1. An aircraft airspeed angle of attack estimation method based on augmented generalized Kalman filtering is characterized in that: The following steps are involved: S1: Obtain gyroscope angular velocity, three-axis acceleration, and altimeter data, project the three-axis acceleration components in the aircraft coordinate system, and then perform back-projection correction. Calculate an initial airspeed estimate based on the altitude difference, calculate an initial angle of attack estimate based on the three-axis acceleration components and the angle between them and gravity, and generate gyroscope state parameter values. S2: Based on the gyroscope state parameter values, extract the airspeed estimate and angle of attack estimate within the estimation period, and the point-by-point difference between the airspeed estimate and angle of attack estimate within the current period, normalize them separately, and then combine them to construct an input matrix to generate a standardized error matrix value; S3: Calling the standardized error matrix value, extracting time series row vectors in sequence, comparing the slope differences between adjacent time series row vectors, dividing the scoring level intervals, assigning scores to the time series row vectors, extracting the score change trend, selecting the slope increase value corresponding to the score change interval as a compensation factor, adjusting the airspeed and angle of attack in the gyroscope state parameter values, and generating an offset dynamic prediction result; S4: Based on the offset dynamic prediction result, a state variable set is formed, the pitch angular velocity, the rudder deflection amount, and the lift and drag coefficient bias value are used as additional quantities, the slope of each state element is extracted and then calculated with the period interval ratio to generate a structural state change coefficient; S5: calling the airspeed change value and the angle of attack change value in the structural state change coefficient, combining the current observation projection residual and the predicted value, extracting the relevant elements of the state error covariance matrix, setting the weighted ratio, and fusing the data from each source proportionally to generate a joint airspeed angle of attack output; The steps for obtaining the combined airspeed angle of attack output are specifically as follows: S501: calling the airspeed change value in the structural state change coefficient, extracting the transverse covariance element in the state error covariance matrix, calculating the ratio of the observation projection residual to the predicted value, setting the weighted proportional coefficient of the transverse covariance element and the longitudinal covariance element, and generating the covariance matrix weighted proportional coefficient; S502: The weighted proportional coefficient of the covariance matrix is called, and the sum of the absolute differences between the observed projection residual and the predicted value is combined to calculate the square root of the square of the airspeed change increment and the angle of attack change increment using the formula: ； Calculate the weighted fusion factor, integrate the vertical elements of the covariance matrix, and generate the multi-source data fusion weight; in, represents the m-th dimension weighted fusion factor, represents the m-th dimension covariance weighted proportional coefficient, represents the qth group of observation projection residuals, represents the predicted value of group q, represents the time series airspeed change increment, represents the increment of the angle of attack change in the spatial sequence, represents the horizontal elements in the covariance matrix, represents the vertical element in the covariance matrix, is the denominator protection constant; S503: Calling the airspeed change value and the angle of attack change value, calculating the weighted summation result according to the multi-source data fusion weight, performing error threshold boundary truncation processing on the fusion result, and generating a joint airspeed angle of attack output.

2. The method for estimating the airspeed angle of attack of an aircraft based on augmented generalized Kalman filtering according to claim 1, wherein: The standardized error matrix value includes a normalized airspeed error vector, a normalized angle of attack error vector, and an estimation period offset amplitude coefficient. The offset dynamic prediction result is specifically a corrected airspeed prediction value, a corrected angle of attack prediction value, and an airspeed and angle of attack change trend fitting factor. The structural state change coefficient includes an airspeed change slope coefficient, an angle of attack change slope coefficient, and a state response time ratio. The joint airspeed and angle of attack output includes an airspeed estimation fusion result, an angle of attack estimation fusion result, and a residual weighted fusion coefficient.

3. The method for estimating the airspeed angle of attack of an aircraft based on augmented generalized Kalman filtering according to claim 2, wherein: The steps for obtaining the standardized error matrix value are specifically as follows: S201: Based on the gyroscope state parameter value, obtain the airspeed estimate and the angle of attack estimate within the estimation period, calculate the point-by-point difference between the current period airspeed value and the estimated period airspeed estimate, and combine the point-by-point difference between the current period angle of attack value and the estimated period angle of attack estimate to generate a period difference sequence matrix; S202: calling the periodic difference sequence matrix, dividing the multi-periodic point difference values by the corresponding maximum difference value for normalization, combining the normalized airspeed difference values and the angle of attack difference values to construct a matrix, and obtaining a normalized combined input matrix; S203: Based on the normalized combined input matrix, the normalized airspeed difference and angle of attack difference of multiple period points are called, and the difference with the root mean square value and the average value of the corresponding channel is calculated using the formula: ； The standardized error values of multiple period points are obtained by operation and combined to generate the standardized error matrix value; in, , represents the normalized airspeed difference at the qth period point, represents the normalized angle of attack difference at the qth period point, represents the RMS value of the angle of attack difference sequence of the p-th channel in the current estimation cycle, Represents the average value of the normalized combined input matrix of the p-th channel in the current estimation cycle.

4. The method for estimating the airspeed angle of attack of an aircraft based on augmented generalized Kalman filtering according to claim 3, wherein: The steps for obtaining the offset dynamic prediction result are specifically as follows: S301: calling the standardized error matrix value, extracting adjacent time series row vectors according to the row vector index order, calculating the slope difference of corresponding elements between adjacent time series row vectors, taking the absolute value and then the average, and generating a slope difference value; S302: Based on the slope difference values, divide the scoring level intervals according to a preset threshold range, compare each slope difference value with the boundary values of multiple intervals step by step, store the scoring values corresponding to the matching intervals into an array, and generate a scoring vector; S303: Call the rating vector to calculate the incremental difference between adjacent rating values using the formula: ； The compensation factor is obtained by calculation and linearly superimposed with the initial airspeed nominal value and the angle of attack baseline parameter to generate the offset dynamic prediction result; in, represents the compensation factor of the qth time series row vector at time s, Represents the slope difference value of adjacent time series row vectors, represents the rating vector, Represents the difference between the mth adjacent rating changes, represents the length of the rating vector, represents the time interval between adjacent time series row vectors, represents the angle of attack adjustment in the previous iteration cycle, Represents the upper threshold of angle of attack.

5. The method for estimating the airspeed angle of attack of an aircraft based on augmented generalized Kalman filtering according to claim 4, wherein: The steps for obtaining the structural state variation coefficient are specifically as follows: S401: Based on the offset dynamic prediction result, the pitch angular velocity, the control surface deflection, and the lift and drag coefficient offset are called, and the three are aligned with the lateral displacement data and longitudinal acceleration data of the gyroscope state parameters in time series. The three are combined into a multidimensional set of dynamic parameters and basic state parameters to generate an extended set of state variables; S402: extracting the time-domain variation slope of the pitch angular velocity, the time-domain variation slope of the rudder deflection, and the time-domain variation slope of the lift and drag coefficient bias value from the state variable extended set, performing linear fitting on the time-domain variation slopes of the three parameters using the least squares method, and calculating the absolute value of the variation rate to construct a state element slope set; S403: Calling the time-domain variation slope of the pitch angular velocity, the time-domain variation slope of the rudder deflection, and the time-domain variation slope of the lift and drag coefficient bias value in the state element slope set, and calculating the ratio of the absolute value of the multiple slopes to the period interval parameter using the formula: ； By calculating the dynamic trend principal component and the rudder coupling compensation term, the structural state variation coefficient parameters are generated by superposition; in, represents the structural state variation coefficient parameter, Represents the slope value of the i-th element in the state element slope set, represents the period interval parameter, Represents the minimum parameter of the lift-drag coefficient bias value period, Represents the maximum value parameter of the rudder deflection, represents the dynamic coupling factor parameter of the rudder surface, represents the periodic average value parameter of the lift-drag coefficient bias value, Represents the lift-drag stability compensation parameter.

6. The method for estimating the airspeed angle of attack of an aircraft based on augmented generalized Kalman filtering according to claim 1, wherein: The steps for obtaining the gyroscope state parameter value are specifically as follows: S101: Based on the gyroscope angular velocity, three-axis acceleration, and altimeter data, the linear acceleration is obtained and projected into the aircraft coordinate system. Combined with the attitude angle, it is back-projected into the inertial system. The velocity change is calculated based on the altitude difference to generate an initial airspeed estimate. S102: Call the initial airspeed estimate and the longitudinal acceleration component, calculate the angle with the gravity direction, and determine the incident direction and attitude offset using the formula: ； The calculation obtains the estimated value of the longitudinal incident angle, obtains the deviation trend based on the difference between it and the attitude angle, obtains the initial estimate of the angle of attack, and generates the gyroscope state parameter value; in, represents the initial estimate of the angle of attack, Represents the acceleration component in the Z-axis direction in the body coordinate system, represents the acceleration due to gravity, represents the pitch angle, represents the initial estimate of airspeed, Respectively represent the acceleration components of the X-axis and Y-axis directions in the body coordinate system, Represents the altitude values obtained by the altimeter at the current moment and the previous moment respectively.

Images