Aircraft airspeed attack angle estimation method based on augmented generalized Kalman filtering

Abstract

The invention relates to the technical field of aircraft control, in particular to an aircraft airspeed attack angle estimation method based on augmented generalized Kalman filtering, which comprises the following steps: acquiring sensor data and correcting acceleration, calculating initial estimation of an airspeed attack angle, constructing an error matrix and scoring, extracting variation trend adjustment parameters, and calculating a state variation coefficient. And multi-source information is fused to output the combined airspeed angle of attack quantity. According to the method, through acceleration component body coordinate system projection and back projection correction, an airspeed initial estimation value is constructed in combination with a height difference value, a physical matching degree is enhanced, a slope score and trend change extraction are performed to establish a dynamic compensation mechanism, and the adaptive capacity to a non-stationary flight state is improved; the introduction of the pitch angle speed, the control surface deflection and the aerodynamic coefficient bias makes the state evolution more structural, a change coefficient is established through ratio calculation, the coupling of state prediction and maneuvering dynamic is enhanced, the weighted fusion of an observation residual error and a predicted value is realized, and the convergence capability of an estimated value to complex aerodynamic interference is improved.

Description

Technical Field

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

[0002] The technical field of aircraft control includes the improvement of the stability and handling performance of aircraft, mainly involving the monitoring, estimation, and control of important parameters such as the airspeed, attitude, heading, and angle of attack of the aircraft. The core content of aircraft control is to precisely control the flight state of the aircraft through various sensors, measurement devices, and control systems to ensure the safe flight of the aircraft in complex environments. This technical field is widely applied to the design and operation of various aircraft such as civil aviation, military aircraft, and unmanned aerial vehicles, and involves multiple sub-fields such as autopilot, flight control systems, and flight state estimation. The key objective of aircraft control technology is to enhance the autonomy, reliability, and safety of aircraft.

[0003] Among them, the aircraft airspeed and angle of attack estimation method refers to using the sensor information and control algorithms of the aircraft to estimate the airspeed and angle of attack of the aircraft in real time. This technology specifically solves the accuracy problem of measuring the airspeed and angle of attack parameters during the flight of the aircraft, especially in situations where direct measurement is difficult. By processing the data from multiple sensors of the aircraft and combining mathematical models, the estimation method dynamically obtains the airspeed and angle of attack of the aircraft, thereby providing real-time feedback for the automatic control system of the aircraft. These methods generally achieve precise estimation of airspeed and angle of attack by fusing and processing sensor data and using technical means such as Kalman filter and least squares method.

[0004] Relying on the method of static sensor fusion and model estimation is difficult to describe the dynamic nonlinear response during high-speed flight, and the airspeed and angle of attack estimation lags behind the response to sudden attitude changes. The changes in time-domain data are not fully utilized, resulting in the lack of the ability of the estimation result to identify and adjust the trend changes. The structural state modeling is insufficient, and the prediction is inaccurate under multi-axis linkage control or complex aerodynamic environments, affecting the consistency between the estimation and the actual flight state, resulting in deviation of the control response and reduced flight stability. Summary of the Invention

[0005] The purpose of the present invention is to solve the deficiencies existing in the prior art, and propose an aircraft airspeed and angle of attack estimation method based on augmented extended Kalman filter.

[0006] To achieve the above objective, the present invention adopts the following technical solution: An aircraft airspeed and angle of attack estimation method based on augmented extended Kalman filter, including the following steps:

[0007] S2: Based on the gyroscope state parameter values, extract the airspeed and angle of attack estimations within the estimation period, calculate the point-by-point differences with the current period values, normalize them respectively, and then combine them to construct an input matrix, generating a standardized error matrix value;

[0008] S3: Invoke the standardized error matrix value, sequentially extract the time-series row vectors, compare the slope differences between adjacent vectors, set the scoring level division intervals, assign scores to each vector, extract the scoring change trend, and select the slope increase value corresponding to the scoring change interval as the compensation factor to adjust the airspeed and angle of attack in the initial state parameter values, generating a deviation dynamic prediction result;

[0009] S4: Based on the deviation dynamic prediction result, form a state variable set, invoke the pitch angular velocity, control surface deflection amount, and lift-drag coefficient offset values as additional quantities, extract the slopes of each state element, and calculate the ratio with the period interval to generate a structural state change coefficient;

[0010] S5: Invoke the airspeed and angle of attack change values in the structural state change coefficient, combine the current observed projection residual and predicted value, extract the relevant elements of the state error covariance matrix to set the weighting ratio, and fuse the data from each source according to the ratio to generate a combined airspeed and angle of attack output quantity.

[0011] 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 deviation 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 combined airspeed and angle of attack output quantity includes an airspeed estimation fusion result, an angle of attack estimation fusion result, and a residual weighted fusion coefficient.

[0012] As a further solution of the present invention, the specific steps for obtaining the standardized error matrix value are as follows:

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

[0014] S202: Invoke the period difference sequence matrix, respectively normalize the multi-period point differences by dividing them by the corresponding maximum differences, combine the normalized airspeed differences and angle of attack differences to construct a matrix, and obtain a normalized combined input matrix;

[0015] S203: According to the normalized combined input matrix, invoke the multi-period point normalized airspeed differences and angle of attack differences, calculate the differences from the root mean square value and average value of the corresponding channels, and use the formula:

[0016] ;

[0017] Operate to obtain the standardized error value of the multi-cycle points and combine to generate the standardized error matrix value;

[0018] Among them, , represents the normalized airspeed difference of the q-th cycle point, represents the normalized angle of attack difference of the q-th cycle point, represents the root mean square 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.

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

[0020] 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 the corresponding elements between adjacent vectors, take the absolute value and then take the average value to generate a slope difference value;

[0021] S302: Based on the slope difference value, divide the scoring grade interval according to the preset threshold range, compare each difference value with the multi-interval boundary values step by step, and store the scoring value corresponding to the matching interval into an array to generate a scoring vector;

[0022] S303: Call the scoring vector, calculate the incremental difference of adjacent scoring values, and use the formula:

[0023] ;

[0024] Operate to obtain the compensation factor value, linearly superimpose the compensation factor with the initial airspeed nominal value and the angle of attack baseline parameter to generate the offset dynamic prediction result;

[0025] Among them, represents the compensation factor of the q-th time-series row vector at time s, represents the slope difference value of adjacent time-series row vectors, represents the scoring vector, represents the m-th adjacent scoring change difference, represents the length of the scoring vector, represents the time interval of adjacent time-series row vectors, represents the angle of attack adjustment amount of the previous iteration cycle, represents the upper limit of the angle of attack threshold.

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

[0027] S401: Based on the lateral dynamic prediction result, call the pitch angular velocity parameter, the rudder surface deflection parameter, and the lift-drag coefficient offset value parameter, align the three with the lateral displacement data and longitudinal acceleration data of the original state parameters in time series, and merge them into a multi-dimensional set of dynamic parameters and basic state parameters to generate an extended set of state variables;

[0028] S402: Extract the time-domain change slope of the pitch angular velocity parameter, the time-domain change slope of the rudder surface deflection parameter, and the time-domain change slope of the lift-drag coefficient offset value parameter from the extended set of state variables, perform linear fitting on the time series curves of the three parameters by the least squares method, and calculate the absolute value of the change rate to construct a set of state element slopes;

[0029] S403: Call the pitch angular velocity time-domain change slope, the rudder surface deflection time-domain change slope, and the lift-drag coefficient offset value time-domain change slope in the set of state element slopes, calculate the ratio of the multi-slope absolute value to the period interval parameter, and use the formula:

[0030] ;

[0031] Generate the structural state change coefficient parameter by calculating the dynamic trend principal component and the rudder surface coupling compensation term and superimposing them;

[0032] wherein, represents the structural state change coefficient parameter, represents the slope value of the i-th element in the set of state element slopes, represents the period interval parameter, represents the minimum value parameter within the period of the lift-drag coefficient offset value, represents the maximum value parameter of the rudder surface deflection, represents the rudder surface dynamic coupling factor parameter, represents the average value parameter of the period of the lift-drag coefficient offset value, represents the lift-drag stability compensation amount parameter.

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

[0034] S501: Call the airspeed change value in the structural state change coefficient, extract the lateral 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 lateral covariance element and the longitudinal covariance element, and generate the covariance matrix weighted proportional coefficient;

[0035] S502: Call the covariance matrix weighted proportionality coefficient, and combine the sum of the absolute differences between the observed projection residuals and the predicted values to calculate the square root of the sum of the squares of the airspeed change increment and the angle of attack change increment, using the formula:

[0036] ;

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

[0038] Among them, represents the m-th dimensional weighted fusion factor, represents the m-th dimensional covariance weighted proportionality coefficient, represents the q-th group of observed projection residuals, represents the q-th group of predicted values, represents the time series airspeed change increment, represents the spatial series angle of attack change increment, represents the transverse element in the covariance matrix, represents the longitudinal element in the covariance matrix, is the denominator protection constant;

[0039] 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, and perform error threshold boundary truncation processing on the fusion result to generate the combined airspeed angle of attack output.

[0040] As a further solution of the present invention, the method further includes:

[0041] S1: Obtain the gyroscope angular velocity, three-axis acceleration, and altimeter data, perform back-projection correction on the components of the acceleration projected in the body coordinate system, calculate the initial airspeed estimate in combination with the height difference, calculate the initial angle of attack estimate based on the angle between the acceleration component and gravity, and generate the gyroscope state parameter value;

[0042] The gyroscope state parameter value is specifically the body attitude angular rate, acceleration correction residual, and initial angle of attack estimate.

[0043] As a further solution of the present invention, the acquisition steps of the gyroscope state parameter value are specifically:

[0044] S101: Based on the gyroscope angular velocity, three-axis acceleration, and altimeter data, obtain the linear acceleration, project it onto the body coordinate system, back-project it to the inertial system in combination with the attitude angle, and calculate the speed change based on the height difference to generate the initial airspeed estimate;

[0045] S102: Call the initial airspeed estimate and the longitudinal acceleration component, calculate the angle with the gravity direction, and judge the incident direction and attitude deviation, using the formula:

[0046] ;

[0047] Calculate the estimated value of the longitudinal incident angle through operations, obtain the offset trend based on the difference from the attitude angle, and get the initial estimated value of the angle of attack;

[0048] Among them, represents the initial estimated value of the angle of attack, represents the acceleration component in the Z-axis direction in the body coordinate system, represents the gravitational acceleration, represents the pitch angle, represents the initial estimated value of the airspeed, respectively represent the acceleration components in the X-axis and Y-axis directions in the body coordinate system, respectively represent the altitude values obtained by the altimeter at the current moment and the previous moment;

[0049] S103: Invoke the differences in the initial estimated value of the angle of attack and the changes in pitch and roll angular velocities, combine with the lateral acceleration to judge the disturbance change amount, and generate the gyroscope state parameter value.

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

[0051] In the present invention, by performing projection and back-projection corrections of the acceleration components in the body coordinate system, and constructing the initial estimated value of the airspeed in combination with the altitude difference, the physical matching degree of the original estimation is enhanced. The normalization processing of the period difference makes the data input more consistent. The slope scoring and trend change extraction establish a dynamic compensation mechanism, improving the adaptability to non-steady flight states. The introduction of pitch angular velocity, control surface deflection, and aerodynamic coefficient offset makes the state evolution more structured. By calculating the ratio, a change coefficient is established to strengthen the coupling between state prediction and maneuver dynamics. The weighted fusion of the observation residual and the predicted value improves the convergence ability of the estimation to complex aerodynamic disturbances. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0053] Figure 2 is a flowchart of the steps for obtaining the gyroscope state parameter value of the present invention;

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

[0055] Figure 4 is a flowchart of the steps for obtaining the offset dynamic prediction result of the present invention;

[0056] Figure 5 is a flowchart of the steps for obtaining the structure state change coefficient of the present invention;

[0057] Figure 6 This is a flowchart for obtaining the combined airspeed and angle of attack output of the present invention. Specific embodiments

[0058] In order to make the objectives, technical solutions and advantages of the present invention clearer, 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 used to limit the present invention.

[0059] In the description of the present invention, it should be understood that the orientation or positional relationship indicated by the terms "length", "width", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", etc. is based on the orientation or positional relationship shown in the accompanying drawings. These are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and thus should not be construed as limiting the present invention. In addition, in the description of the present invention, "a plurality of" means two or more unless otherwise specifically defined.

[0060] Embodiment 1

[0061] Please refer to Figure 1 , the present invention provides a technical solution: an aircraft airspeed and angle of attack estimation method based on augmented extended Kalman filter, including the following steps:

[0062] S1: Obtain gyroscope angular velocity, three-axis acceleration and altimeter data, perform back-projection correction on the projection of the acceleration components in the body coordinate system, calculate the initial airspeed estimate in combination with the height difference, calculate the initial angle of attack estimate based on the angle between the acceleration component and gravity, and generate gyroscope state parameter values;

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

[0064] S3: Call the standardized error matrix value, sequentially extract the time-series row vectors, compare the slope differences between adjacent vectors, set the scoring level division interval, assign scores to each vector and then extract the scoring change trend, select the slope increase value corresponding to the scoring change interval as the compensation factor, and adjust the airspeed and angle of attack in the initial state parameter values to generate a deviation dynamic prediction result;

[0065] S4: Based on the deviation dynamic prediction result, form a state variable set, call the pitch angular velocity, control surface deflection amount and lift-drag coefficient offset value as additional quantities, extract the slope of each state element and calculate it with the ratio of the period interval to generate a structural state change coefficient;

[0066] S5: Call the airspeed and angle of attack change values in the structural state change coefficient, combine the current observed projection residual and the predicted value, extract the relevant elements of the state error covariance matrix to set the weighting ratio, and fuse the data from each source according to the ratio to generate the combined airspeed and 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 estimation value. The standardized error matrix values include the normalized airspeed error vector, normalized angle of attack error vector, and estimation period offset amplitude coefficient. The offset dynamic prediction results are specifically the corrected airspeed prediction value, corrected angle of attack prediction value, and airspeed and angle of attack change trend fitting factor. The structural state change coefficient includes the airspeed change slope coefficient, angle of attack change slope coefficient, and state response time ratio. The combined airspeed and angle of attack output includes the airspeed estimation fusion result, angle of attack estimation fusion result, and residual weighting fusion coefficient.

[0068] Please refer to Figure 2 , and the specific steps for obtaining the gyroscope state parameter values are as follows:

[0069] S101: Based on the gyroscope angular velocity, three-axis acceleration, and altimeter data, obtain the linear acceleration, project it onto the body coordinate system, and back-project it to the inertial system in combination with the attitude angle. Calculate the speed change based on the height difference to generate the initial airspeed estimate;

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

[0071] 1. Remove the gravity component: Through the attitude angles , (obtained by integrating the gyroscope and fusing with the Kalman filter), calculate the components of gravity on each axis , , to obtain the linear acceleration , , ;

[0072] 2. Project to the inertial system: Through the rotation matrix (calculated from the attitude angle), project to the inertial system acceleration , where ;

[0073] 3. Integrate the speed change: Combine the height difference Calculate the vertical speed (time interval: 1 s). Initial estimated airspeed (corrected to 25 m / s after smoothing filter).

[0074] S102: Call the initial estimated airspeed and the longitudinal acceleration component, calculate the angle with the gravity direction, judge the incident direction and attitude offset, and use the formula:

[0075] ;

[0076] Obtain the estimated value of the longitudinal incident angle through the operation, obtain the offset trend based on the difference with the attitude angle, and get the initial estimated angle of attack;

[0077] Among them, represents the initial estimated angle of attack, represents the acceleration component in the Z-axis direction in the body coordinate system, represents the gravitational acceleration, represents the pitch angle, represents the initial estimated airspeed, respectively represent the acceleration components in the X-axis and Y-axis directions in the body coordinate system, respectively represent the altitude values obtained by the altimeter at the current moment and the previous moment;

[0078] In the process of calculating the initial estimated angle of attack, first call the initial estimated airspeed , collect the acceleration component in the Z-axis of the body coordinate system (measured by a triaxial accelerometer, the original data is , after removing the gravity component ), pitch angle (calculated by integrating the angular velocity of the gyroscope, the angular velocity , the integration time , obtained after error compensation), gravitational acceleration , the acceleration components in the X-axis and Y-axis are respectively , (obtained after coordinate projection), the current altitude , the previous altitude (measured by a barometric altimeter, the sampling interval is ).

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

[0080] 1. Calculate ;

[0081] 2. Calculate ;

[0082] 3. Divide by Obtained 。

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

[0084] 1. Calculate ;

[0085] 2. Calculate 。

[0086] After adding the two terms and multiplying by :

[0087] 1. Calculate ;

[0088] 2. Multiply by Obtained ;

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

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

[0091] Table 1: Table of parameters related to the initial airspeed estimate

[0092] Parameter Value Unit #timg# 25 m / s #timg# 12.5 m / s² #timg# 2 m

[0093] As shown in Table 1, the initial airspeed estimate and related parameters are obtained through sensor data and coordinate transformation. The formula quantifies the initial angle of attack estimate through the ratio relationship between longitudinal acceleration and lateral acceleration, combined with the rate of change of height. The result shows that the estimated angle of attack value exceeds the range of conventional flight states (typical value <15°), and the offset trend needs to be corrected by combining the attitude angle difference.

[0094] S103: Call the differences in the initial angle of attack and the changes in pitch and roll angular velocities, and combine the lateral acceleration to judge the disturbance change amount, and generate the gyroscope state parameter value.

[0095] Call the initial angle of attack estimate , the pitch angular velocity , the roll angular velocity (after differential calculation), the lateral acceleration , and perform the following operations:

[0096] 1. Calculate the angular velocity difference ;

[0097] 2. Judge the disturbance threshold: If and (the threshold is set based on the statistical data of the aircraft's stable state), then it is determined that there is a disturbance;

[0098] 3. Generate gyroscope state parameters: disturbance intensity , state parameters (for control law correction).

[0099] The advantage of the formula is that by fusing the longitudinal acceleration component and the lateral acceleration ratio relationship, and dynamically adjusting the angle-of-attack estimation weight in combination with the height change rate, the problem that single-sensor data is vulnerable to interference is solved.

[0100] Please refer to Figure 3 , and the specific steps for obtaining the standardized error matrix value are as follows:

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

[0102] Based on the gyroscope state parameter values , set the estimation period (including 5 data points, sampling interval 1 s), call the current-period airspeed value (measured by the pitot tube and smoothed by Kalman filter), the airspeed estimation value within the estimation period (weighted prediction through state parameters and historical data), the current-period angle-of-attack value , the angle-of-attack estimation value within the estimation period , and perform the following operations:

[0103] 1. Calculate the point-by-point difference of airspeed ( to 5), to obtain ;

[0104] 2. Calculate the point-by-point difference of angle of attack , to obtain ;

[0105] 3. Arrange the difference sequence in matrix form according to the periodic points:

[0106] ;

[0107] Table 2: Gyroscope state parameters and associated data table

[0108] Parameter Value Unit #timg# 1.638 None #timg# 24.8-25.0 m / s #timg# 31.2-31.7 Degree

[0109] 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 the same as the number of periodic points.

[0110] S202: Call the matrix of cycle difference sequences, and normalize the multi-cycle point differences by dividing them by their corresponding maximum differences respectively. Combine the normalized airspeed differences and angle of attack differences to construct a matrix, and obtain the normalized combined input matrix;

[0111] Call the matrix of cycle difference sequences and extract the column of airspeed differences , the column of angle of attack differences , and perform the following operations:

[0112] 1. Calculate the maximum value of airspeed differences , the maximum value of angle of attack differences ;

[0113] 2. Normalize the airspeed differences: , and obtain ;

[0114] 3. Normalize the angle of attack differences: , and obtain ;

[0115] 4. Combine them into a normalized input matrix:

[0116] ;

[0117] S203: According to the normalized combined input matrix, call the normalized airspeed differences and angle of attack differences of multi-cycle points, and calculate the differences from the root mean square values and average values of the corresponding channels, using the formula:

[0118] ;

[0119] Obtain the standardized error values of multi-cycle points through operations, and combine them to generate the standardized error matrix values;

[0120] Among them, .

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

[0122] 1. For channel 1 (airspeed differences):

[0123] The root mean square value ;

[0124] The average value ;

[0125] 2. For channel 2 (angle of attack differences):

[0126] The root mean square value ;

[0127] Average value ;

[0128] Substitute into the formula to calculate the standardized error value (taking as an example):

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

[0130] 2. ;

[0131] 3. That is ;

[0132] Similarly, calculate the error values of all periodic points to generate a standardized error matrix:

[0133] ;

[0134] The advantage of the formula is that by introducing the root mean square value and average value of the channels, the weight of the normalized difference is dynamically adjusted to eliminate the interference of parameters with different dimensions on error evaluation. The result shows that the standardized error matrix can quantify the fluctuation characteristics of airspeed and angle of attack estimation within a period, providing an input basis for subsequent control compensation.

[0135] Please refer to Figure 4 , and the specific steps for obtaining the offset dynamic prediction result are as follows:

[0136] 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 the corresponding elements between adjacent vectors, take the absolute value and then take the average to generate a slope difference value;

[0137] Call the standardized error matrix value, set the adjacent time-series row vector index to 4 (the number of rows of the matrix is 5, and the adjacent rows are separated by 1 second), and extract the error matrix as follows:

[0138] ;

[0139] Perform the following operations:

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

[0141] For and , the slope difference of the first channel (airspeed error) is , the slope difference of the second channel (angle of attack error) is , and the average value after taking the absolute value is ;

[0142] For and , the slope difference of the first channel is , the second channel is , and the mean value is ;

[0143] For and , the slope difference of the first channel is , the second channel is , and the mean value is ;

[0144] For and , the slope difference of the first channel is , the second channel is , and the mean value is ;

[0145] 2. Generate a sequence of slope difference values: .

[0146] Table 3: Example Table of Standardized Error Matrix

[0147] Timing row index Airspeed error channel Angle of attack error channel 1 0.3 0.4 2 0.2 0.3 3 0.1 0.5 4 0.4 0.2 5 0.3 0.6

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

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

[0150] Based on the sequence of slope difference values , the preset scoring grade intervals are:

[0151] Low fluctuation (0 ≤ difference value < 0.1): Score 1;

[0152] Medium fluctuation (0.1 ≤ difference value < 0.2): Score 2;

[0153] High fluctuation (difference value ≥ 0.2): Score 3;

[0154] Perform step-by-step comparison:

[0155] 1. The difference value of 0.1 belongs to the medium fluctuation interval, score 2;

[0156] 2. The difference value of 0.15 belongs to the medium fluctuation interval, score 2;

[0157] 3. The difference value of 0.3 belongs to the high fluctuation interval, score 3;

[0158] 4. The difference value of 0.25 belongs to the high fluctuation range, and the score is 3;

[0159] Generate a score vector: 。

[0160] Table 4: Matching Table of Score Level Intervals and Difference Values

[0161] Difference value range Score Instance difference value Matching result [0, 0.1) 1 0.05 Score 1 [0.1, 0.2) 2 0.15 Score 2 [0.2, ∞) 3 0.25 Score 3

[0162] As shown in Table 4, the score level is achieved by gradually comparing with the interval boundary values, and the difference value is mapped one-to-one with the score result.

[0163] S303: Call the score vector, calculate the incremental difference between adjacent score values, using the formula:

[0164] ;

[0165] Operate to obtain the compensation factor value, linearly superimpose the compensation factor with the initial airspeed nominal value and the angle of attack baseline parameter to generate a deviation dynamic prediction result;

[0166] Among them, represents the compensation factor of the q-th time series row vector at time s, represents the slope difference value between adjacent time series row vectors, represents the score vector, represents the m-th adjacent score change difference, represents the length of the score vector, represents the time interval between adjacent time series row vectors, represents the angle of attack adjustment amount in the previous iteration cycle, represents the upper limit of the angle of attack threshold.

[0167] Call the score vector , the angle of attack adjustment amount in the previous iteration cycle, the upper limit of the angle of attack threshold (set according to the aerodynamic performance of the aircraft), the time interval , calculate the adjacent score increment difference :

[0168] 1. ;

[0169] 2. ;

[0170] 3. ;

[0171] Substitute into the formula to calculate the compensation factor (taking as an example):

[0172] ;

[0173] Similarly, calculate the compensation factors for other cycle points and generate the results: , and combine the compensation factor with the initial airspeed nominal value , angle of attack baseline for linear superposition (weight 0.1) to obtain the predicted offset:

[0174] Airspeed offset: ;

[0175] Angle of attack offset: .

[0176] The advantage of the formula is that by integrating the scoring increment and the angle of attack historical adjustment amount, the adaptability of the compensation factor to time sensitivity and aerodynamic constraints is dynamically corrected. The result shows that the compensation factor is adjusted according to the ratio of the scoring level and the angle of attack threshold, effectively suppressing the overshoot risk.

[0177] Please refer to Figure 5 , and the specific steps for obtaining the structural state change coefficient are as follows:

[0178] S401: Based on the dynamic prediction result of the offset, call the pitch angular velocity parameter, the rudder surface deflection parameter, and the lift-drag coefficient offset value parameter, align the three with the lateral displacement data and longitudinal acceleration data of the original state parameters in time series, and merge them into a multi-dimensional set of dynamic parameters and basic state parameters to generate an extended set of state variables;

[0179] Based on the dynamic prediction result of the offset (airspeed offset , angle of attack offset ), call the pitch angular velocity parameter (measured by the gyroscope), the rudder surface deflection parameter (fed back by the actuator), the lift-drag coefficient offset value (calibrated through wind tunnel experiments), the lateral displacement data (measured by GPS), the longitudinal acceleration data (collected by the accelerometer), align them in time series (sampling interval ), and merge them into a multi-dimensional set:

[0180] ;

[0181] Table 5: Example table of the extended set of multi-dimensional state variables

[0182] Pitch angular velocity (rad / s) Rudder deflection (°) Lift-drag offset value Lateral displacement (m) Longitudinal acceleration (m / s²) 0.10 2.1 0.05 0.20 2.1 0.15 2.3 -0.02 0.25 2.3 0.12 2.5 0.03 0.30 2.0 0.18 2.0 -0.01 0.18 2.4

[0183] As shown in Table 5, the multi-dimensional set is aligned through timestamps, integrating the dynamic prediction parameters and the basic state parameters.

[0184] S402: Extract the time-domain change slope of the pitch angular velocity parameter, the time-domain change slope of the rudder surface deflection amount parameter, and the time-domain change slope of the lift-drag coefficient offset value parameter from the extended set of state variables. Perform linear fitting on the time-series curves of the three parameters by the least squares method, and calculate the absolute value of the change rate to construct a set of state element slopes;

[0185] Extract the pitch angular velocity, rudder surface deflection amount, and lift-drag offset value parameters from the extended set, and perform least squares linear fitting on the time-series data of each parameter (4 cycle points):

[0186] 1. Pitch angular velocity slope:

[0187] Data points: , ;

[0188] Slope calculation: ;

[0189] 2. Rudder surface deflection amount slope:

[0190] Data points: ;

[0191] Slope calculation: ;

[0192] 3. Lift-drag offset value slope:

[0193] Data points: ;

[0194] Slope calculation: ;

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

[0196] S403: Call the time-domain change slope of the pitch angular velocity, the time-domain change slope of the rudder surface deflection amount, and the time-domain change slope of the lift-drag coefficient offset value in the set of state element slopes, calculate the ratio of the multi-slope absolute value to the period interval parameter, and use the formula:

[0197] ;

[0198] Generate the structural state change coefficient parameter by calculating the dynamic trend principal component and the rudder surface coupling compensation term;

[0199] Among them, represents the structural state change coefficient parameter, represents the slope value of the i-th element in the set of state element slopes, represents the period interval parameter, represents the minimum value parameter within the lift-drag coefficient offset value period, Parameter representing the maximum value of rudder surface deflection, Parameter representing the dynamic coupling factor of the rudder surface, Parameter representing the periodic average value of the lift and drag coefficient offset, Parameter representing the lift and drag stability compensation amount.

[0200] Call the slope set , the period interval parameter , the minimum value within the lift and drag coefficient offset period (extracted from the fourth row of Table 1), the maximum value of rudder surface deflection , the dynamic coupling factor of the rudder surface (calibrated according to the rudder effectiveness curve), the periodic average value of the lift and drag coefficient , the stability compensation amount (to prevent the denominator from being zero), substitute into the formula:

[0201] ;

[0202] Calculate step by step:

[0203] 1. Calculate the sum of the squares of the slopes: ;

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

[0205] 3. The denominator term: ;

[0206] 4. The result of the first term: ;

[0207] 5. The numerator of the second term: ;

[0208] 6. The denominator of the second term: ;

[0209] 7. The result of the second term: ;

[0210] 8. The final .

[0211] The advantage of the formula is that by integrating the dynamic trend principal component and the rudder surface coupling effect, it quantifies the comprehensive influence of aerodynamics and control. The result shows that the structural state change coefficient is significantly dominated by the rudder surface deflection, and the rudder effectiveness compensation needs to be adjusted preferentially.

[0212] Please refer to Figure 6 , the specific steps for obtaining the combined airspeed and angle of attack output are as follows:

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

[0214] Call the structural state change coefficient in the airspeed change value , extract the lateral covariance element in the state error covariance matrix (corresponding to the lateral displacement error), the longitudinal covariance element (corresponding to the longitudinal acceleration error), the observed projection residual (the difference between the GPS measurement and the prediction), the predicted value , and perform the following operations:

[0215] 1. Calculate the ratio of the observed residual to the predicted value: ;

[0216] 2. Set the weighted proportionality coefficient of the lateral and longitudinal covariances ;

[0217] 3. Generate the set of covariance matrix weighted proportionality coefficients: .

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

[0219] Covariance type Value Unit Lateral displacement error 0.12 m2 Longitudinal acceleration error 0.08 m2 / s2

[0220] As shown in Table 6, the weighted proportionality coefficient is directly calculated through the ratio of the covariance elements.

[0221] S502: Call the covariance matrix weighted proportionality coefficient, combine the total absolute difference between the observed projection residual and the predicted value, calculate the square root of the sum of the squares of the airspeed change increment and the angle of attack change increment, and use the formula:

[0222] ;

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

[0224] Among them, represents the m-th dimensional weighted fusion factor, represents the m-th dimensional covariance weighted proportionality coefficient, represents the q-th group of observed projection residuals, represents the q-th group of predicted values, represents the time series airspeed change increment, represents the spatial series angle of attack change increment, represents the lateral element in the covariance matrix, Represents the longitudinal elements in the covariance matrix, is the denominator protection constant;

[0225] Call the weighted proportionality coefficient , the sum of the absolute differences between the observed residuals and the predicted values , the increment of airspeed change , the increment of angle of attack change (converted to radians: ), the denominator protection constant , substitute into the formula to calculate the weighted fusion factor (taking as an example):

[0226] ;

[0227] Similarly calculate , , generate the multi-source data fusion weights: .

[0228] 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 weights, and perform error threshold boundary truncation processing on the fusion result to generate the combined airspeed and angle of attack output.

[0229] Call the airspeed change value , the angle of attack change value , the fusion weight , perform weighted summation:

[0230] ;

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

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

[0233] If it exceeds, truncate it to the boundary value.

[0234] The benefit of the formula is that by fusing covariance ratios and dynamic increments, it balances the reliability of horizontal and vertical data. The result shows that the weighted factor is dominated by the covariance ratio, ensuring the attenuation of weights in high-error dimensions.

[0235] The above are only the preferred embodiments of the present invention, and do not limit the present invention in other forms. Any person skilled in the art may use the technical content disclosed above to make changes or modifications into equivalent embodiments with equivalent changes and apply them to other fields. However, as long as it does not depart from the technical solution content of the present invention, any simple modification, equivalent change, and modification made to the above embodiments based on the technical essence of the present invention still fall within the protection scope of the technical solution of the present invention.

Claims

1. An aircraft airspeed and angle of attack estimation method based on augmented extended Kalman filter, characterized in that It includes the following steps: S2: Based on the gyroscope state parameter values, extract the airspeed and angle of attack estimations within the estimation period, calculate the point-by-point differences with the current period values, normalize them respectively, and then combine them to construct an input matrix, and generate a standardized error matrix value; S3: Call the standardized error matrix value, sequentially extract the time-series row vectors, compare the slope differences between adjacent vectors, set the scoring level division intervals, assign scores to each vector, then extract the scoring change trend, and select the slope increase value corresponding to the scoring change interval as the compensation factor to adjust the airspeed and angle of attack in the initial state parameter values to generate a deviation dynamic prediction result; S4: Based on the deviation dynamic prediction result, form a state variable set, call the pitch angular velocity, rudder deflection amount, and lift and drag coefficient offset values as additional quantities, calculate the ratio of the slope of each state element to the period interval to generate a structural state change coefficient; S5: Call the airspeed and angle of attack change values in the structural state change coefficient, combine the current observed projection residual and prediction value, extract the relevant elements of the state error covariance matrix to set the weighting ratio, and fuse the data from each source according to the ratio to generate a combined airspeed and angle of attack output.

2. The method for estimating the airspeed and angle of attack of an aircraft based on the augmented extended Kalman filter 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 deviation 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 combined 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 and angle of attack of an aircraft based on augmented extended Kalman filter according to claim 2, wherein The specific steps for obtaining the standardized error matrix value are as follows: S201: Based on the gyroscope state parameter values, obtain the airspeed estimation and angle of attack estimation within the estimation period, calculate the point-by-point differences between the current period airspeed value and the airspeed estimation in the estimation period, and combine the point-by-point differences between the current period angle of attack value and the angle of attack estimation in the estimation period to generate a period difference sequence matrix; S202: Call the period difference sequence matrix, respectively normalize the multi-period point differences by dividing them by the corresponding maximum difference, combine the normalized airspeed differences and angle of attack differences to construct a matrix, and obtain a normalized combined input matrix; S203: According to the normalized combined input matrix, call the multi-period point normalized airspeed differences and angle of attack differences, and calculate the differences from the root mean square value and average value of the corresponding channels. Use the formula: ； Perform operations to obtain the standardized error values of multi-period points, and combine them to generate a standardized error matrix value; Among them, , represents the normalized airspeed difference of the q-th cycle point, represents the normalized angle of attack difference of the q-th cycle point, represents the root mean square 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 and angle of attack of an aircraft based on the augmented extended Kalman filter according to claim 3, wherein The specific steps for obtaining the deviation dynamic prediction result are as follows: S301: Call the standardized error matrix value, extract adjacent time-series row vectors according to the row vector index order, calculate the slope differences of the corresponding elements between adjacent vectors, take the absolute value and then take the average to generate a slope difference value; S302: Based on the slope difference value, divide the scoring level intervals according to the preset threshold range, compare each difference value with the multi-interval boundary values step by step, and store the scoring values corresponding to the matching intervals in an array to generate a scoring vector; S303: Call the scoring vector, calculate the incremental difference between adjacent scoring values, using the formula: ； Calculate to obtain the compensation factor value, linearly superimpose the compensation factor with the initial airspeed nominal value and the angle of attack baseline parameter to generate a deviation dynamic prediction result; Among them, represents the compensation factor of the q-th timing row vector at time s, represents the slope difference value between adjacent timing row vectors, represents the scoring vector, represents the m-th adjacent scoring change difference, represents the length of the scoring vector, represents the time interval between adjacent timing row vectors, represents the angle of attack adjustment amount in the previous iteration cycle, represents the upper limit of the angle of attack threshold.

5. The method for estimating the airspeed and angle of attack of an aircraft based on augmented extended Kalman filter according to claim 4, characterized in that The specific steps for obtaining the structural state change coefficient are as follows: S401: Based on the deviation dynamic prediction result, call the pitch angular velocity parameter, the rudder deflection amount parameter, and the lift-drag coefficient offset value parameter, align the three with the lateral displacement data and longitudinal acceleration data of the original state parameter in time series, and merge them into a multi-dimensional set of dynamic parameters and basic state parameters to generate a state variable expansion set; S402: Extract the time-domain change slope of the pitch angular velocity parameter, the time-domain change slope of the rudder deflection amount parameter, and the time-domain change slope of the lift-drag coefficient offset value parameter from the state variable expansion set, linearly fit the time series curves of the three parameters by the least squares method, and calculate the absolute value of the change rate to construct a state element slope set; S403: Call the pitch angular velocity time-domain change slope, the rudder deflection amount time-domain change slope, and the lift-drag coefficient offset value time-domain change slope in the state element slope set, calculate the ratio of the multi-slope absolute value to the period interval parameter, and use the formula: ； Generate the structural state change coefficient parameter by calculating the dynamic trend principal component and the rudder coupling compensation term and superimposing them; Among them, represents the structural state change coefficient parameter, represents the slope value of the i-th element in the set of state element slopes, represents the period interval parameter, represents the minimum value parameter within the period of the lift-drag coefficient offset value, represents the maximum value parameter of the rudder surface deflection amount, represents the rudder surface dynamic coupling factor parameter, represents the period average value parameter of the lift-drag coefficient offset value, represents the lift-drag stability compensation amount parameter.

6. The method for estimating the airspeed and angle of attack of an aircraft based on augmented extended Kalman filter according to claim 5, wherein The specific steps for obtaining the combined airspeed and angle of attack output are as follows: S501: Call the airspeed change value in the structural state change coefficient, extract the lateral covariance element in the state error covariance matrix, calculate the ratio of the observation projection residual to the predicted value, set the weighted proportionality coefficient of the lateral covariance element and the longitudinal covariance element to generate the covariance matrix weighted proportionality coefficient; S502: Call the covariance matrix weighted proportionality coefficient, combine the sum of the absolute differences between the observation projection residual and the predicted value, calculate the square root of the sum of the squares of the airspeed change increment and the angle of attack change increment, and use the formula: ； Calculate the weighted fusion factor, integrate the longitudinal elements of the covariance matrix to generate the multi-source data fusion weight; Among them, represents the m-th dimensional weighted fusion factor, represents the m-th dimensional covariance weighted proportionality coefficient, represents the q-th group of observed projection residuals, represents the q-th group of predicted values, represents the time series airspeed change increment, represents the spatial series angle of attack change increment, represents the lateral element in the covariance matrix, represents the longitudinal element in the covariance matrix, is the denominator protection constant; 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, and perform error threshold boundary truncation processing on the fusion result to generate the combined airspeed and angle of attack output.

7. The method for estimating the airspeed and angle of attack of an aircraft based on the augmented extended Kalman filter according to claim 6, wherein The method further includes: S1: Obtain the gyroscope angular velocity, three-axis acceleration, and altimeter data, perform back-projection correction on the projection of the acceleration components in the body coordinate system, calculate the initial airspeed estimate in combination with the height difference, calculate the initial angle of attack estimate based on the angle between the acceleration component and gravity, and generate the gyroscope state parameter value; The gyroscope state parameter value is specifically the body attitude angle rate, the acceleration correction residual, and the initial angle of attack estimate value.

8. The method for estimating the airspeed and angle of attack of an aircraft based on the augmented extended Kalman filter according to claim 7, characterized in that, The specific steps for obtaining the gyroscope state parameter value are as follows: S101: Based on the gyroscope angular velocity, three-axis acceleration, and altimeter data, obtain the linear acceleration, project it onto the body coordinate system, back-project it to the inertial system in combination with the attitude angle, and calculate the speed change based on the height difference to generate the initial airspeed estimate; S102: Call the initial airspeed estimate and the longitudinal acceleration component, calculate the angle with the gravity direction, and judge the incident direction and attitude deviation, using the formula: ； Calculate to obtain the estimated value of the longitudinal incident angle, obtain the offset trend according to the difference from the attitude angle, and obtain the initial estimated value of the angle of attack; Among them, represents the initial estimated value of the angle of attack, represents the acceleration component in the Z-axis direction in the body coordinate system, represents the gravitational acceleration, represents the pitch angle, represents the initial estimated value of the airspeed, respectively represent the acceleration components in the X-axis and Y-axis directions in the body coordinate system, respectively represent the altitude values obtained by the altimeter at the current moment and the previous moment; S103: Call the difference in the change of the initial estimated value of the angle of attack and the pitch and roll angular velocities, combine with the lateral acceleration to judge the disturbance change amount, and generate the gyroscope state parameter value.

Images