Open AccessArticle

Optimization-Assisted Filter for Flow Angle Estimation of SUAV Without Adequate Measurement

Ziyi Wang

Jie Li

¹,

Chang Liu

^2,*,

Yu Yang

²,

Juan Li

¹,

Xueyong Wu

¹,

Yachao Yang

² and

Bobo Ye

School of Mechatronical Engineering, Beijing Institute of Technology, Beijing 100081, China

Yangtze Delta Region Academy of Beijing Institute of Technology, JIAXING, No. 1940, East North Road, Youchegang Town, Xiuzhou District, Jiaxing 314019, China

Author to whom correspondence should be addressed.

Drones 2024, 8(12), 758; https://doi.org/10.3390/drones8120758

Submission received: 9 November 2024 / Revised: 9 December 2024 / Accepted: 10 December 2024 / Published: 15 December 2024

Download

Browse Figures

Versions Notes

Abstract

The accurate estimation of flow angles is crucial for enhancing flight performance and aircraft safety. Flow angles of fixed-wing small unmanned aerial vehicles (SUAVs) are more vulnerable due to their low airspeed. Current flow angle measurement devices have not been widely implemented in SUAVs due to their substantial cost and size constraints. Moreover, there are no general estimation methods suitable for SUAVs based on their rudimentary sensor suite. This study presents a generalized optimization-assisted filter estimation (OAFE) method for estimating the relative velocity and flow angles of fixed-wing SUAVs based on a standard sensor suite. This OAFE method mainly consists of a cubature Kalman filter and an optimizer. The filter serves as the main loop with which to generate flow angles in real time by fusing the acceleration, angular rate, attitude, and airspeed. Without flow angle measurements, the optimizer generates approximate aerodynamic derivatives, which serve as pseudo-measurements with which to refine the performance of the filter. The results demonstrate that the estimated angle of attack and side slip angle displayed root mean square errors of around 0.11° and 0.24° in the simulation. The feasibility was also verified in field tests. The OAFE method does not require flow angle measurements, the prior acquisition of aerodynamic parameters, or model training, making it suitable for quick deployment on different SUAVs.

Keywords:

angle of attack; flow angle; estimation; fixed-wing SUAV; cubature Kalman filter

1. Introduction

In recent years, fixed-wing small unmanned aerial vehicles (SUAVs) have evolved rapidly, and are frequently used due to their affordability and convenience. The angle of attack (AoA or

α

) and side slip angle (SSA or

β

)—the two flow angles—are important parameters that describe the air data of fixed-wing SUAVs. The accurate estimation of flow angles is crucial for enhancing the flight performance and safety of an aircraft. Moreover, advanced flight controllers [1], system identification [2], and trajectory prediction [3] usually require mapping between the flow angles and the aerodynamic force.

Compared with general aviation, the flow angles of SUAVs are more vulnerable to wind field and maneuvering due to their smaller size, lighter weight, and lower airspeed. However, obtaining flow angles for SUAVs is even more challenging due to the lack of economical, compact measurement devices.

Extensive research has been conducted on flow angle acquisition for SUAVs. Currently, flow angle acquisition methods for SUAVs can be classified into the following three primary categories: (1) flow angle measurement techniques; (2) kinematic-based estimation methods; and (3) model-aided estimation methods.

Flow angle measurement techniques mainly involve the use of multi-hole probes, vane probes, multi-pitot tubes, and distributed pressure sensors [4]. The main feature of this technique is the direct acquisition of flow angles based on sufficient sensor information, generally without the use of inertial measurements. Only a few studies have fused inertial information in order to reduce the noise of vane probe measurements [5]. On the one hand, calibration constitutes a primary research focus concerning multi-hole probes and vane probes [6]. A question arises about the accuracy of such probes at low speed, and their calibration must be performed in a wind tunnel. Furthermore, probes with mobile vanes are too large and delicate to be easily integrated. On the other hand, the main purpose of research on multi-pitot tubes and distributed pressure sensors is mostly to reduce the costs and embed the sensors within the surfaces of UAVs [7]. The mathematical models of these devices are intricate. As a result, multiple machine learning algorithms have been utilized for modeling these flow angle sensors in recent studies [7,8]. This data-driven sensor modeling approach requires a large dataset. Liu et al. developed dimensionless input and output neural networks to improve the generalization in the entire flight envelope and effectively reduce the need for training data [9]. Nevertheless, due to differences in aerodynamic profiles, measurement techniques using distributed sensors require dedicated training data for specific aircraft. The result is that the equipment, especially for SUAVs, is rarely available at a low cost [4].

Kinematic-based estimation methods rely primarily on the wind triangle, incorporating the ground speed, wind field, and attitude of SUAVs [10,11]. SUAVs equipped with a standard sensor suite provide access to ground speed, attitude, and 1-D airspeed data. However, these sensor data are insufficient to support the simultaneous estimation of wind field and flow angles using recursive filters. The challenge of insufficient system observability was identified by Wenz et al. [12], who employed an optimization estimation method in their subsequent work. However, this also causes issues with real-time performance [13]. Therefore, many kinematic-based methods still adopt the 2D horizontal wind assumption [14] and the zero-sideslip assumption, which reduce the dimensionality of variants [15]. Some studies in general aviation assume that there is no wind [16]. However, assumptions of the kinematic-based method are not well suited for low-speed SUAVs, as they lead to larger flow angle errors.

The model-aided method does not require information related to wind field and flow angle measurements but does require accurate aerodynamic parameters of the aircraft. These parameters can be either classical derivatives or a neural network [17], as long as the relationship between air data and aerodynamic forces is mapped. The model-aided method that first appeared in [18] estimated the AoA based on aerodynamic stability derivatives using only data from accelerometers and control surface positions. Tian et al. realized real-time estimation of flow angles based on the aerodynamic derivatives of SUAVs; the estimation accuracy of the AoA reached a standard deviation of 0.96° [5,19,20]. If air data are accessible, a custom neural network model can be developed using a data-driven approach to represent the input–output relationship between the measurements and air data. After the custom neural network is obtained, it can be used to replace the air data systems [21,22]. However, a neural network model trained on flight data is interdependent with model-aided flow angle estimation. As a result, this method was used for synthetic airspeed, AoA, and SSA estimation for analytical redundancy when considering the presence of noise in the measurement devices [23]. Finally, model-aided methods can only be used for a specific UAV if the aerodynamic parameters of that particular UAV are available. The model-aided method is not readily applicable to all SUAVs, as obtaining the aerodynamic parameters of a specific SUAV is a complex task.

Overall, all three categories of methods have significant limitations when applied to SUAVs. The method proposed in this study draws on the model-aided method, which utilizes the aerodynamics of SUAVs. The main difference is that the proposed method does not require prior knowledge of the specific parameters of the aerodynamic model. The main contributions of this study are summarized as follows:

(1) An optimization-assisted filter estimation (OAFE) method is proposed for estimating flow angles in fixed-wing SUAVs. The OAFE relies solely on a standard sensor suite comprising a 1D pitot tube, GPS, and an IMU, without dependence on the information of wind field or aerodynamic parameters of a specific SUAV. Therefore, it can be quickly deployed to different SUAVs;

(2) A cubature Kalman Filter is designed as the main loop to ensure real-time performance. “2D horizontal wind” and “zero-sideslip” assumptions are not adopted in the filter, which makes the method applicable to SUAVs at low speed. An optimizer loop was developed to provide real-time aerodynamic derivatives as pseudo-measurements. Although these real-time derivatives are noisy, this is to be expected since there is no flow angle information during the optimization process;

(3) To evaluate the OAFE, both 3D flight simulations and field tests were conducted. The simulation utilized an open-source flight dynamics model and a professional autopilot. The effectiveness of the estimation of relative velocity, flow angles, aerodynamic derivatives, and specific force is validated in this study.

The organization of this paper is as follows. Section 2 presents the mathematical models of flow angle estimation and introduces the theoretical frameworks employed in this paper. The proposed methods are discussed in detail in Section 3. In Section 4 and Section 5, the proposed method is evaluated through simulation and field tests conducted on a fixed-wing SUAV. Section 6 summarizes this study and indicates future work.

2. Problem Formulation

2.1. Kinematics Model

The definition of flow angles is illustrated in Figure 1, where

α

denotes the angle of attack and

β

denotes the side slip angle. Generally, the relative velocity (also referred to as the velocity relative to air) of the aircraft to the surrounding airmass is labeled as

v_{r}

. The relative velocity is labeled as

v_{r}^{b} = {[u, v, w]}^{T}

in the body frame, as well as

v_{r}^{n}

in the navigation frame. Below, the superscript

b

denotes values in the body frame,

n

denotes the navigation frame, and

a

denotes the wind frame.

Once

v_{r}^{b}

is known, the airspeed

V_{a}

α

β

, and dynamic pressure

\bar{q}

can be calculated by Equation (1) as follows:

\begin{array}{l} V_{a} = ∥ [u, v, w] ∥ \\ α = \tan^{- 1} (w / u) \\ β = \sin^{- 1} (v / V_{a}) \\ \bar{q} = 0.5 ρ V_{a}^{2} \end{array}

(1)

where

ρ

is the air density.

Wind significantly affects the air data of SUAVs flying at low speed. The ground velocity

v_{g}

is not a good approximation of

v_{r}

[24]. As shown in Equation (2), the so-called “wind triangle” denotes the relationship between

v_{g}

, wind velocity over ground

v_{w}

, and

v_{r}

in the navigation frame [25]. To convert this to the body frame, the rotation matrix

C_{n}^{b}

(refer to Appendix A.1) is introduced, as follows:

\begin{array}{l} v_{r}^{n} = v_{g}^{n} - v_{w}^{n} \\ v_{r}^{b} = C_{n}^{b} (v_{g}^{n} - v_{w}^{n}) \end{array}

(2)

Equation (2) provides a direct way to compute

v_{r}^{b}

. However, the standard sensor suite for SUAVs does not provide access to the real-time 3D wind field. Available engineering methods make zero-sideslip and horizontal wind field assumptions [26], which are practical approaches for macroscopic wind field estimation, but not beneficial for the precise estimation of

v_{r}^{b}

2.2. Dynamics Model

The SUAV’s aerodynamic force

F_{a}

is given in [27], which, when decomposed in the wind frame, is denoted as follows:

F_{a} = {[- D, Y, - L]}^{T} = \bar{q} \bar{S} [\begin{matrix} - (C_{D_{0}} + C_{D_{α}} α^{2}) \\ C_{Y_{β}} β + {C_{Y}}_{δ_{r}} \\ - (C_{L_{0}} + C_{L_{α}} α + C_{L_{q}} q \bar{c} / 2 U_{0} + C_{L_{δ_{e}}} δ_{e}) \end{matrix}]

(3)

where

\bar{S}

is the wing area,

C_{*_{*}}

are aerodynamic derivatives,

\bar{c}

is the mean aerodynamic chord,

U_{0}

is the trim airspeed,

q

is the pitch angular rate, and

δ_{e}

and

δ_{r}

are the elevator and rudder deflections.

For low-speed SUAVs, neglecting the effects of the Earth’s rotation and revolution, the dynamics are as follows:

m {\dot{v}}_{g}^{b} = m (g^{b} + f^{b}) = G^{b} + F_{T}^{b} + C_{a}^{b} F_{a}

(4)

where

m

is the mass,

g^{b}

is the gravitational acceleration,

f^{b}

is the specific force,

G^{b}

is gravity,

F_{T}^{b}

is the thrust, and

C_{a}^{b}

is the rotation matrix from the wind frame to the body frame given in Appendix A.1.

Considering SUAVs equipped with an IMU,

f^{b}

can be measured directly. After cancelling out gravity, Equation (4) is represented as follows:

m f^{b} = F_{T}^{b} + C_{a}^{b} F_{a}

(5)

Equations (3) and (5) show the functional relation between aerodynamic derivatives and flow angles. Flow angles can be estimated on the premise that the exact aerodynamic derivatives of the SUAV are known [28].

2.3. Dynamics Simplification

SUAVs typically use non-vectorized propulsion. When decomposed in the body frame,

F_{T}^{b}

has no component along the

y

and

z

axes. For small flow angle flights, Equations (3) and (5) yield the relationship among

[C_{Y_{β}}, {C_{Y}}_{δ_{r}}, C_{L_{0}}, C_{L_{α}}, C_{L_{q}}, C_{L_{δ_{e}}}]

v_{r}^{b}

, and

f^{b}

as follows:

\begin{array}{l} f^{b, y} = \bar{q} \bar{S} (C_{Y_{β}} β + {C_{Y}}_{δ_{r}}) \cos β / m \\ f^{b, z} = - \bar{q} \bar{S} (C_{L_{0}} + C_{L_{α}} α + C_{L_{q}} \frac{\bar{c}}{2 U_{0}} q + C_{L_{δ_{e}}} δ_{e}) \cos α / m \end{array}

(6)

where

f^{b, y}

and

f^{b, z}

denote components of

f^{b}

along the

y

and

z

-axes.

2.4. Cubature Kalman Filter (CKF)

CKF is an accurate nonlinear filter that can be applied to solve high-dimensional nonlinear filtering problems with minimal computational efforts [29]. It uses a series of cubature points to propagate a priori and a posteriori statistical characteristic. As a result, the calculation of Jacobian matrix is avoided. In order to accommodate the fast maneuvers of aircraft, we used the square-root version of the CKF for its numerical stability. The square-root CKF algorithm is summarized as follows [30]:

(1) Time update

The predicted state

{\hat{x}}_{k |k - 1}

and the square-root factor of the predicted error covariance

S_{k |k - 1}

, are computed using the estimated state

{\hat{x}}_{k - 1 |k - 1}

, the square root of state error covariance

S_{k - 1 |k - 1}

, and input

u_{k - 1}

, as follows:

X_{i, k - 1 |k - 1} = S_{k - 1 |k - 1} λ_{i} + {\hat{x}}_{k - 1 |k - 1}

(7)

X_{i, k |k - 1}^{*} = f (X_{i, k - 1 |k - 1}, u_{k - 1})

(8)

{\hat{x}}_{k |k - 1} = \frac{1}{2 n} \sum_{i = 1}^{2 n} X_{i, k |k - 1}^{*}

(9)

χ_{k |k - 1}^{*} = \frac{1}{\sqrt{2 n}} [\begin{matrix} X_{1, k |k - 1}^{*} - {\hat{x}}_{k |k - 1} & \dots & X_{2 n, k |k - 1}^{*} - {\hat{x}}_{k |k - 1} \end{matrix}]

(10)

S_{k |k - 1} = Tria ([\begin{matrix} χ_{k |k - 1}^{*} & S_{Q, k - 1} \end{matrix}])

(11)

where

n

is the dimension of the state

x

i = 1,2, \dots, 2 n

λ_{i} = \sqrt{n} {[1]}_{i}

and

f (\cdot)

is the process equations.

T r i a (\cdot)

denotes a general triangularization algorithm, which returns a lower triangular matrix.

(2) Measurement update

The updated state

{\hat{x}}_{k |k}

and

S_{k | k}

are estimated using

{\hat{x}}_{k |k - 1}

S_{k |k - 1}

, measurement

z_{k}

, and the square root of the measurement error covariance

S_{R, k}

, as follows:

X_{i, k |k - 1} = S_{k |k - 1} λ_{i} + {\hat{x}}_{k |k - 1}

(12)

Z_{i, k |k - 1} = h (X_{i, k |k - 1}, u_{k})

(13)

{\hat{z}}_{k |k - 1} = \frac{1}{2 n} \sum_{i = 1}^{2 n} Z_{i, k |k - 1}

(14)

Z_{k | k - 1} = \frac{1}{\sqrt{2 n}} [\begin{matrix} Z_{1, k |k - 1} - {\hat{z}}_{k |k - 1} & \dots & Z_{2 n, k |k - 1} - {\hat{z}}_{k |k - 1} \end{matrix}]

(15)

S_{z z, k |k - 1} = Tria ([\begin{matrix} Z_{k | k - 1} & S_{R, k} \end{matrix}])

(16)

χ_{k |k - 1} = \frac{1}{\sqrt{2 n}} [\begin{matrix} X_{i, k |k - 1} - {\hat{x}}_{k |k - 1} & \dots & X_{2 n, k |k - 1} - {\hat{x}}_{k |k - 1} \end{matrix}]

(17)

P_{x z, k |k - 1} = χ_{k |k - 1} Z_{k |k - 1}^{T}

(18)

W_{k} = (P_{x z, k |k - 1} / S_{z z, k |k - 1}^{T}) / S_{z z, k |k - 1}

(19)

{\hat{x}}_{k |k} = {\hat{x}}_{k |k - 1} + W_{k} (z_{k} - {\hat{z}}_{k |k - 1})

(20)

S_{k ∣ k} = Tria ([\begin{matrix} χ_{k ∣ k - 1} - W_{k} Z_{k ∣ k - 1} & W_{k} S_{R, k} \end{matrix}])

(21)

where

h (\cdot)

is the measurement equation.

3. OAFE Method

This study proposes an OAFE method for flow angle estimation comprising a filter loop and an optimizer loop (shown in Figure 2). The OAFE combines the benefits of both the filter and optimizer, with the filter serving as the main loop to generate flow angles.

A CKF filter was employed with low computational complexity and high frequency to capture states (

x_{k}

) changes caused by the SUAV maneuvering. The filter incorporates

{\tilde{z}}_{p, k}

from the optimizer loop as pseudo-measurements, thereby reducing cumulative errors and ensuring real-time outputs. The subscript

k

in Figure 2 denotes values at time

t_{k}

Angular rates

ω

, specific force

f

, pitot airspeed

u_{p}

, Euler angles

ξ

, and

[{\tilde{C}}_{Y_{β}} {\tilde{C}}_{L}_{α}]

are assumed to be available. Where

ω

f

u_{p}

ξ

are bias-corrected and provided by the integrated navigation system,

[{\tilde{C}}_{Y_{β}} {\tilde{C}}_{L}_{α}]

are real-time derivatives identified by the optimizer loop.

The optimizer loop generates aerodynamic derivatives at a frequency

f_{e}

, based on a sequence (SEQ) of past sensor data prior to the moment

t_{k}

(as shown in Algorithm 1). The SEQ of the data is collected within a moving horizon, the length of which is

n_{o} + 1

Algorithm 1 Optimizer for generating pseudo-measurement

3.1. Filter Loop

In this study, SUAVs are treated as nonlinear systems, as follows:

\begin{matrix} \begin{array}{l} x_{k} = f (x_{k - 1}, u_{k - 1}) + ν_{k - 1} \\ z_{k} = h (x_{k}, u_{k},) + w_{k} \end{array} \end{matrix}

(22)

where

T

is the sample time interval,

v_{k - 1}

and

w_{k}

are the uncorrelated process and measurement Gaussian noise with zero means.

In this study, the filter loop is designed based on the dynamics and kinematics of the SUAV. Based on the CKF introduced in Section 2.4, the state (output)

x

, input

u

, and measurement

z

vectors are given as follows:

x = {[u, v, w, {C_{L}}_{0}, {C_{L}}_{α}, C_{L_{q}}, C_{L_{δ_{e}}}, C_{Y_{β}}]}^{T}

(23)

u = {[p, q, r, ϕ, θ, ψ, f^{b, x}, δ_{e}, δ_{r}]}^{T}

(24)

\begin{matrix} z = {[u_{p}, f^{b, y}, f^{b, z}, {\tilde{C}}_{L}_{α}, {\tilde{C}}_{Y_{β}}]}^{T} \end{matrix}

(25)

Section 3.1.1 and Section 3.1.2 are explicit descriptions of the designed filter in the context of this study.

3.1.1. State Dynamic and Measurement

(1) State Dynamic Equations

The aerodynamic derivatives in state

x

use a static propagation model, as follows:

[\dot{{C_{L}}_{0}}, \dot{{C_{L}}_{α}}, \dot{C_{L_{q}}}, \dot{C_{L_{δ_{e}}}}, \dot{C_{Y_{β}}}] = 0

(26)

The dynamic equation of the

v_{r}^{b}

in state

x

can be obtained by differentiating Equation (2). Due to the high filter frequency, the wind field varies very little in a single period [31]; therefore,

{\dot{v}}_{r}^{b}

can be represented as follows:

\begin{array}{l} \dot{v_{r}^{b}} = \dot{C_{n}^{b}} [v_{g}^{n} - v_{w}^{n}] + C_{n}^{b} [\dot{v_{g}^{n}} - {\dot{v}}_{w}^{n}] \\ \approx C_{n}^{b} (ω_{b \leftarrow n}^{n} \times) v_{r}^{n} + C_{n}^{b} \dot{v_{g}^{n}} \\ = - C_{n}^{b} (ω_{n \leftarrow b}^{n} \times) C_{b}^{n} C_{n}^{b} v_{r}^{n} + \dot{v_{g}^{b}} \\ = - (ω_{n \leftarrow b}^{b} \times) v_{r}^{b} + \dot{v_{g}^{b}} \end{array}

(27)

where

ω_{b \leftarrow n}

represents the angular rate of the navigation frame relative to the body frame while

ω_{n \leftarrow b}

represents the opposite;

C_{n}^{b}

is an orthogonal matrix, and the derivative of the rotation matrix

{\dot{C}}_{n}^{b}

is given in [32].

Combining Equations (4) and (27), the UAV relative velocity is presented as follows:

\dot{v_{r}^{b}} = - (ω_{i \leftarrow b}^{b} \times) v_{r}^{b} + (G^{b} + F_{T}^{b} + C_{a}^{b} F_{a}) / m

(28)

where

ω_{i \leftarrow b}^{b} \approx ω_{n \leftarrow b}^{b}

denotes the angular rates measured by IMU, the following expression is succinctly expressed as

ω^{b}

(ω^{b} \times)

is the antisymmetric matrix consisting of

[p, q, r]

and is denoted by

(ω^{b} \times) = [\begin{matrix} 0 & r & - q \\ - r & 0 & p \\ q & - p & 0 \end{matrix}]

Thus, the dynamic equations can be represented as follows:

\begin{array}{l} \dot{x} = {[\dot{u}, \dot{v}, \dot{w}, \dot{{C_{L}}_{0}}, \dot{{C_{L}}_{α}}, \dot{C_{L_{q}}}, \dot{C_{L_{δ_{e}}}}, \dot{C_{Y_{β}}}]}^{T} \\ = - [\begin{matrix} - (ω^{b} \times) v_{r}^{b} \\ 0_{(5 \times 1)} \end{matrix}] + [\begin{matrix} {C_{n}^{b} g^{n}}_{(3 \times 1)} \\ 0_{(5 \times 1)} \end{matrix}] + [\begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} f^{b, x} \\ \bar{q} \bar{S} (C_{Y_{β}} β + {C_{Y}}_{δ_{r}} δ_{r}) \cos β / m \\ - \bar{q} \bar{S} (C_{L_{0}} + C_{L_{α}} α + C_{L_{q}} \frac{\bar{c}}{2 V_{a}} q + C_{L_{δ_{e}}} δ_{e}) \cos α / m \end{matrix} \end{matrix} \\ 0_{(5 \times 1)} \end{matrix} \end{matrix}] \end{array}

(29)

where

f^{b, x}

denotes the component of

f^{b}

along the x-axes.

The cubature points propagate in Equation (8) according to the dynamic equations described in Equation (29).

(2) Measurement equations

The functional relation between

[u_{p}, {\tilde{C}}_{Y_{β}}, {\tilde{C}}_{L}_{α}]

in corrector

z

and the state vector

x

is evident. Combined with Equation (6), the measurement equation

h (\cdot)

in Equations (13) and (22) is precisely defined by Equation (30). where

δ_{r}

q

, and

δ_{e}

are disregarded as they change at high frequency.

h (x, u) = [\begin{matrix} \begin{matrix} \begin{matrix} u \\ \bar{q} \bar{S} C_{Y_{β}} β \cos β / m \\ - \bar{q} \bar{S} (C_{L_{α}} α + C_{L_{0}}) \cos α / m \end{matrix} \\ C_{Y_{β}} \end{matrix} \\ C_{L_{α}} \end{matrix}]

(30)

3.1.2. Filtering Noise

The IMU data used for the filters are bias-corrected. The estimation of these biases is conducted in the integrated navigation module, in which the bias of

ω^{b}

and

f^{b}

are generally modeled as a first-order Markov process [15]. The main consideration in this study is the IMU noise, which affects both the state time update and the measurement update.

(1) Process error covariance

S_{Q, k - 1}

in Equation (11) is a square-root factor of process error covariance

Q_{k - 1}

, such that

Q_{k - 1} = S_{Q, k - 1} S_{Q, k - 1}^{T}

Q_{k - 1}

is the covariance of process noise

ν_{k - 1}

in Equation (22).

The process error covariance of

v_{r}^{b}

is denoted as

Q_{v, k - 1}

, and that of

[{C_{L}}_{0}, \dots, C_{Y_{β}}]

is denoted as

Q_{C, k - 1}

, such that

Q_{k - 1}

is expressed as follows:

Q_{k - 1} = [\begin{matrix} Q_{v, k - 1 (3 \times 3)} & 0_{(3 \times 5)} \\ 0_{(5 \times 3)} & Q_{C, k - 1 (5 \times 5)} \end{matrix}]

(31)

Errors in the time update of

v_{r}^{b}

are mainly due to input

u

, arising from two terms (

A

) and (

B

), shown as follows:

v_{r, k}^{b} = v_{r, k - 1}^{b} + \int_{t_{k - 1}}^{t_{k}} [\begin{matrix} \underset{(A)}{\underset{⏟}{- (ω^{b} \times) \cdot v_{r}^{b}}} \\ + \underset{(B)}{\underset{⏟}{C_{n}^{b} g^{n}}} + (F_{a}^{b} + F_{T}^{b}) / m \end{matrix}] d t

(32)

In this study, the noises associated with (

A

) and (

B

) are considered to be independent of each other within a filter period and can be evaluated separately. Then, the process error covariance of

v_{r}^{b}

can be expressed as follows:

Q_{v, k - 1} = Q_{v, k - 1}^{A} + Q_{v, k - 1}^{B}

(33)

The propagation residual

e_{A}

v_{r}^{b}

caused by the term (

A

) can be expressed as follows:

e_{A} = \int_{t_{k - 1}}^{t_{k}} - (δ_{ω^{b}} \times) \cdot v_{r, k - 1}^{b} d t

(34)

where

δ_{ω^{b}}

denotes the noise vector of the gyroscope.

According to the derivation in Appendix A.2,

Q_{v, k - 1}^{A}

can be expressed as follows:

Q_{v, k - 1}^{A} = (v_{r, k - 1}^{b} \times) σ_{ω^{b}}^{2} {(v_{r, k - 1}^{b} \times)}^{T} \cdot T

(35)

where

σ_{ω^{b}}^{2}

represents the angular rate variance, which is detailed in Section 4.

For small and low-speed SUAVs, the variation of

g^{n}

can be ignored and

Q_{v, k - 1}^{B}

mainly results from attitude errors as Equation (36). The derivation is detailed in Appendix A.2, as follows:

Q_{v, k - 1}^{B} = C_{n}^{b} (g^{n} \times) σ_{ξ}^{2} (g^{n} \times) C_{b}^{n} \cdot T

(36)

where

σ_{ξ}^{2}

represents the Euler angles variance, which is detailed in Section 4.

(2) Measurement noise

S_{R, k}

in Equation (16) is a square-root factor of the measurement error covariance

R_{k}

, such that

R_{k} = S_{R, k} S_{R, k}^{T}

R_{k}

is the covariance of measurement noise

r_{k}

, and is expressed as follows:

R_{k} = [\begin{matrix} {R_{s}}_{(3 \times 3)} & 0_{(3 \times 2)} \\ 0_{(2 \times 3)} & {R_{p}}_{(2 \times 2)} \end{matrix}]

(37)

where

R_{s}

is constructed from the sensor noise and

R_{p}

denotes the error covariance of real-time derivatives from the optimizer loop at the time

t_{k}

R_{s}

is represented as follows:

R_{s} = d i a g ([σ_{u_{p}}^{2}, σ_{f^{b, y}}^{2}, σ_{f^{b, z}}^{2}])

(38)

where the symbol “

d i a g (a)

” stands for a diagonal matrix with

a

on the main diagonal;

[σ_{u_{p}}^{2}, σ_{f^{b, y}}^{2}, σ_{f^{b, z}}^{2}]

denote the variance of

[u_{p}, f^{b, y}, f^{b, z}]

, and the associated sensor noise is detailed in Section 4.

As mentioned earlier, derivatives exhibit a high level of noise. Real-time [

{\tilde{C}}_{Y_{β}} {\tilde{C}}_{L}_{α}

] generated from the optimizer loop in flight tests are shown in Figure 3. To use them as pseudo-measurements, their noise level

R_{p}

is detailed in the design of the optimizer loop (shown in Section 3.2.3).

3.2. Optimizer Loop

The main difficulty faced by the optimizer loop is that the flow angles are unknown since they are exactly the states to be estimated by the OAFE.

Based on an SEQ of past sensor data collected over multiple periods, the optimizer loop determines a set of aerodynamic derivatives that satisfy the following two conditions:

(1) The SUAV flight states satisfy the motion-continuity rule (Equation (28)) [33];

(2) The air data parameters (including

v_{r}^{b}

) and aerodynamic parameters determine the SUAV’s forces (Equation (6)) [34].

3.2.1. Objective Function

Assuming there is a functional relation between the input and output data, this optimization aims to determine the functional relation that minimizes the error, as follows:

Y = W (Θ) + w

(39)

where

Y

represents the output and

W (\cdot)

represents the input.

Θ

is the unknown variable to be determined.

More precisely,

Y

and

Θ

are defined as follows:

\{\begin{array}{l} Y = {[u_{p}, f^{b, y}, f^{b, z}]}^{T} \\ Θ = {[\tilde{u}, \tilde{v}, \tilde{w}, {\tilde{C}}_{Y_{β}}, {\tilde{C}}_{L_{0}}, {\tilde{C}}_{L_{α}}]}^{T} \end{array}

(40)

where

[\tilde{u}, \tilde{v}, \tilde{w}]

denote the components of

{\tilde{v}}_{r}^{b}

to be determined and

[{\tilde{C}}_{Y_{β}}, {\tilde{C}}_{L_{0}}, {\tilde{C}}_{L_{α}}]

denote the derivatives to be determined.

W (\cdot),

in Equation (39), can be represented as Equation (41). The aerodynamic forces due to

[δ_{r}, q, δ_{e}]

are also ignored, as follows:

\underset{Y}{\underset{⏟}{[\begin{matrix} u_{p} \\ f^{b, y} \\ f^{b, z} \end{matrix}]}} = \underset{W (Θ)}{\underset{⏟}{[\begin{matrix} [\begin{matrix} 1,0, 0 \end{matrix}] {[\tilde{u}, \tilde{v}, \tilde{w},]}^{T} \\ F (\tilde{u}, \tilde{v}, \tilde{w}, {\tilde{C}}_{Y_{β}}, {\tilde{C}}_{L_{0}}, {\tilde{C}}_{L_{α}}) \\ G (\tilde{u}, \tilde{v}, \tilde{w}, {\tilde{C}}_{Y_{β}}, {\tilde{C}}_{L_{0}}, {\tilde{C}}_{L_{α}}) \end{matrix}]}} + w

(41)

where

F (\cdot)

and

G (\cdot)

denote the functional relation among

\tilde{u}

\tilde{v}

\tilde{w}

{\tilde{C}}_{Y_{β}}

{\tilde{C}}_{L_{0}}

{\tilde{C}}_{L_{α}}

, and

[f^{b, y}, f^{b, z}]

in Equation (6);

w

denotes the associated measurement noise.

At time

t_{k}

, when extending Equation (41) to

n_{o} + 1

periods, there are historical measurement sequences, as follows:

\{\begin{array}{l} u_{p, k - n_{o}}, u_{p, k - n_{o} + 1}, \dots, u_{p, k} \\ f_{k - n_{o}}^{b, y}, f_{k - n_{o} + 1}^{b, y}, \dots, f_{k}^{b, y} \\ f_{k - n_{o}}^{b, z}, f_{k - n_{o} + 1}^{b, z}, \dots, f_{k}^{b, z} \end{array},

and the variable sequences are

Θ_{k - n_{o}}, Θ_{k - n_{o} + 1}, \dots Θ_{k}

At time

t_{k}

, subscript

[k - n_{o}, k - n_{o} + 1, \dots, k]

is marked as

[1, 2, \dots, n_{o} + 1]

for easy expression. Finally, the optimizer loop finds the unknown

Θ_{1}, Θ_{2}, \dots Θ_{n_{o} + 1}

, by minimizing the objective function at time

t_{k}

, as follows:

J (Θ) = \min_{Θ_{i}} \sum_{i = 1}^{n_{o} + 1} K \cdot {[\begin{matrix} Y_{i} - W (Θ_{i}) \\ Θ_{1} - {\hat{Θ}}_{p r e} \end{matrix}]}^{T} [\begin{matrix} Y_{i} - W (Θ_{i}) \\ Θ_{1} - {\hat{Θ}}_{p r e} \end{matrix}]

(42)

where

K

denotes the weight vector; one available

K

for testing in this study is

K = [1, 4, 1, 0, 0, 0, 1.5, 1, 1]

;

{\hat{Θ}}_{p r e}

denotes the previous estimates at time

t_{k - n_{o}}

from the filter loop.

Based on our tests using the selected sensors with a data acquisition frequency of 10 Hz, it was considered that a time window width of

n_{o} = 5

in Equation (42) was appropriate. As the acquisition frequency increases, it is advisable to correspondingly increase the value of

n_{o}

. Similarly, if the sensor accuracy improves,

n_{o}

can be also increased.

In this study, the Levenberg–Marquardt algorithm was employed to conduct the optimization task. To ensure prompt optimization outcomes, a maximum of 30 iterations was set for each optimization period.

3.2.2. Improving Optimization Efficiency

(1) Variable dimensionality reduction

At time

t_{k}

[{\tilde{C}}_{Y_{β}}

{\tilde{C}}_{L_{0}}

{\tilde{C}}_{L_{α}}]

Θ_{i}

across the time window are treated as constants. In contrast,

{\tilde{v}}_{r, i}^{b}

is rapidly changing. In order to reduce the dimensionality, the unknown variables

{\tilde{v}}_{r, i}^{b}

Θ_{i}

are treated as functions of

{\tilde{v}}_{r, 1}^{b}

In addition, since the online resolution of Equation (42) uses an optimization algorithm, the optimization process computes

{\tilde{v}}_{r, i}^{b}

multiple times, utilizing distinct

[{\hat{u}}_{1}, {\tilde{v}}_{1}, {\tilde{w}}_{1}]

. In order to avoid a high number of repetitive recursive operations, the increment of

{\tilde{v}}_{r}^{b}

is used to construct the unknown

{\tilde{v}}_{r, i}^{b}

(shown in Figure 4).

Neglecting the variation of the wind field during the time window, the equation is first established in the navigation frame as follows:

{\tilde{v}}_{r, i}^{n} = {\tilde{v}}_{r, 1}^{n} + ∆ v_{i}^{n}, i = 1, 2, \dots, n_{o} + 1

(43)

where

∆ v_{i}^{n}

denotes the increment of

v^{n}

within the time window

[t_{1}, t_{i}]

The discrete form of pre-integration [35] is utilized to approximate

∆ v_{i}^{n}

as follows:

∆ v_{i}^{n} = \sum_{j = 1}^{i - 1} [C_{b}^{n} (j) f_{j}^{b} + g^{n}] \cdot T, i = 1, 2, \dots, n_{o} + 1

(44)

where the script

j

denotes the values at time

t_{j}

and

C_{b}^{n} (j)

represents the rotation matrix at time

t_{j}

As a result, the constructed sequence

{\tilde{v}}_{r, i}^{n}

naturally satisfies the continuity rule expressed in Equation (28). When decomposed in the body frame, Equation (43) is expressed as follows:

{\tilde{v}}_{r, i}^{b} = C_{n}^{b} (i) [C_{b}^{n} (1) {\tilde{v}}_{r, 1}^{b} + ∆ v_{i}^{n}], i = 1, 2, \dots, n_{o} + 1

(45)

(2) Minimize redundant operations

Combining Equations (44) and (45) yields Equation (46). Based on Equation (46), it is possible to directly construct the unknown variables of sequence

{\tilde{v}}_{r, i}^{b}

instead of the repetitive recursive operations. At time

t_{k}

, only

{\tilde{v}}_{1}

and

{\tilde{w}}_{1}

are unknown variables in sequence

{\tilde{v}}_{r, i}^{b}

, as follows:

\begin{array}{l} {\tilde{v}}_{r, i}^{b} = C_{n}^{b} (i) C_{b}^{n} (1) {\tilde{v}}_{r, 1}^{b} + ∆ v_{r, i}^{b}, i = 1, 2, \dots, n_{o} + 1 \\ ∆ v_{r, i}^{b} = C_{n}^{b} (i) \sum_{j = 1}^{i - 1} [C_{b}^{n} (j) f_{j}^{b} + g^{n}] \cdot T \end{array}

(46)

According to Equation (46), the sequence

{\tilde{v}}_{r, i}^{b}

are computed using

C_{n}^{b} (i)

and

∆ v_{r, i}^{b}

. Once

C_{n}^{b} (i)

and

∆ v_{r, i}^{b}

are obtained, they remain constant throughout the optimization process at time

t_{k}

Moreover, Equation (47) is derived from Equation (46). With Equation (47), the sequence

{∆ v_{r, i}^{b} |}_{k + 1}

at time

t_{k + 1}

can be quickly obtained based on the sequence

{∆ v_{r, i}^{b} |}_{k}

at time

t_{k}

. Thus, the optimization task at time

t_{k + 1}

can start quickly, avoiding the repetitive computation of

∆ v_{r, i}^{b}

, as follows:

\begin{matrix} {∆ v_{r, i}^{b} |}_{k + 1} = {C_{n}^{b} (i) |}_{k + 1} \{C_{b}^{n} {(i) |}_{k} \cdot {∆ v_{r, i + 1}^{b} |}_{k} + [{C_{b}^{n} (i - 1) |}_{k + 1} \cdot {f_{i - 1}^{b} |}_{k + 1} - C_{b}^{n} {(1) |}_{k} \cdot {f_{1}^{b} |}_{k}] \cdot T\}, \\ i = 1,2, \dots, n_{o} + 1 \end{matrix}

(47)

when

i = n_{o} + 1

{∆ v_{r, i + 1}^{b} |}_{k} = 0 .

3.2.3. Noises of Selected Derivatives

The sequence

{\tilde{v}}_{r, i}^{b}

can be calculated through Equation (46) at time

t_{k}

, once

[\tilde{v}, \tilde{w}]

are determined. It is easy to calculate sequences

[Q_{1}, \dots, Q_{n_{o} + 1}]

[α_{1}, \dots, α_{n_{o} + 1}]

, and

[β_{1}, \dots, β_{n_{o} + 1}]

from time

t_{k - n}

t_{k}

according to Equation (1).

The error variance is thus expressed as follows:

σ^{2} ({\tilde{C}}_{L}_{α}) = \frac{1}{n_{o} + 1} \sum_{i = 1}^{n_{o} + 1} {((\frac{m f_{i}^{b, z}}{Q_{i} S} - {\tilde{C}}_{L}_{0}) / α_{i} - {\tilde{C}}_{L}_{α})}^{2}

(48)

σ^{2} ({\tilde{C}}_{Y_{β}}) = \frac{1}{n_{o} + 1} \sum_{i = 1}^{n_{o} + 1} {(\frac{m f_{i}^{b, y}}{Q_{i} S β_{i}} - {\tilde{C}}_{L}_{β})}^{2}

(49)

On the other hand, considering that the flight control system counteracts the perturbations required to excite the SUAV dynamics [36], the error covariance of the estimated derivatives is weighted by normalized flow angle variations within the time window. Therefore,

R_{p}

in Equation (37) is set as follows:

R_{p} = [\begin{matrix} σ^{2} ({\tilde{C}}_{L}_{α}) / e^{∆ α} & 0 \\ 0 & σ^{2} ({\tilde{C}}_{L}_{β}) / e^{∆ β} \end{matrix}]

(50)

where

∆ α = \max (α_{k - n}, \dots, α_{k}) - \min (α_{k - n}, \dots, α_{k})

∆ β = \max (β_{k - n}, \dots, β_{k}) - \min (β_{k - n}, \dots, β_{k})

, and

e

is a natural constant.

4. Simulation

In this study, the OAFE method was first tested using a flight simulation. The flight simulation provides a fundamental guarantee for field trials, allows for more scenarios to be simulated, and facilitates more detailed and refined testing of the method in the future. The ideas for evaluations in the simulation are shown in Table 1. Both the flow angles and aerodynamic parameters estimated by the OAFE were evaluated. To enhance readability, the aerodynamic derivatives were synthesized into aerodynamic coefficients for presentation.

4.1. Initialization of OAFE

To start the OAFE, the initial states

x_{0}

, and the initial error covariance

P_{0}

x_{0}

were set as follows:

x_{0} = {[u_{p, 0}, 0, 0, - 2, - 10, 0, 0]}^{T}

(51)

P_{0} = d i a g ([5^{2}, 3^{2}, 3^{2}, 0.01, 0.01, 0.01, 0.01])

(52)

where the

u_{p, 0}

is the pitot airspeed when the filter starts, and

S_{0, 0}

in Section 2.4 is the square root of

P_{0}

Referring to the parameters of onboard sensors, the noise characteristics of the inputs and measurements are shown in Table 2. The values of

ϕ

θ

, and

ψ

provided are reference values; in practice, they are determined by the integrated navigation system and are updated in real time.

The process error covariance of aerodynamic parameters

Q_{C}

in Equation (31) was set as follows:

Q_{C, k - 1} = {0.01}^{2} d i a g ([{0.1}^{2}, 1^{2}, 1^{2}, 1^{2}, 1^{2}])

(53)

4.2. HIL Simulation System

To accurately replicate real-world conditions, a hardware-in-loop (HIL) simulation was implemented. The simulation employed an onboard autopilot (comprising both the hardware Pixhawk and software px4) and an onboard computer system NX3 (running the OAFE) identical to those used in an actual SUAV. A workstation was used for the flight simulation and a control station was used for the mission setup. The hardware configuration is depicted in Figure 5.

We modeled the sensors of the test SUAV in HIL based on the noise characteristics from the real onboard sensors shown in Table 3.

The simulation framework is illustrated in Figure 6. To start the flight simulation, JSBSim first instantiates the SUAV system model and configures the associated parameters. Subsequently, it sets the initial position, attitude, and velocity of the SUAV.

The following processes occurred during each simulation period: ① the JSBSim received commands from the px4 autopilot to control the SUAV’s actuators. Subsequently, it then updated the aerodynamic forces and moments of the model and proceeded with the kinematic solution; ② The px4 autopilot operated on the Pixhawk hardware platform. It acquired sensor data from JSBSim, fused the data for navigation, and generated control commands according to the preset mission; and ③ the OAFE accepted

f^{b}

ξ

ω^{b}

, and

u_{P}

from the px4 and estimated

[u, v, w]

[α, β]

and

[{C_{L}}_{0}, {C_{L}}_{α}, C_{L_{q}}, C_{L_{δ_{e}}}, C_{Y_{β}}]

4.3. Simulation Setup

4.3.1. Simulation Scenes

Wind is an important topic in flight mechanics and the aerodynamic forces of SUAVs are highly sensitive to the wind field. In this study, five different 3D wind fields were set up in the flight simulation. These five wind fields were defined in the navigation frame, including the crosswind, headwind, tailwind, and updraft, as shown in Table 4.

The wind velocity

v_{w}^{n}

was assumed to be the sum of a low-frequency part

v_{s}^{n}

and a high-frequency turbulent part

v_{t}^{n}

. The

v_{s}^{n}

was assumed to be slowly time-varying as

{\dot{v}}_{s}^{n} \approx 0

. The

v_{t}^{n}

is modeled by the Dryden wind model¹⁴ as follows:

v_{t, k + 1}^{n} = v_{t, k}^{n} - T V_{a} (v_{t}^{n} \circ \frac{1}{L_{k}}) + σ_{v_{t}} \circ \sqrt{2 T V_{a} / L_{k}} \circ v_{v_{t}}

(54)

where

L_{k}

is the spatial wavelength,

σ_{v_{t}}

is the noise amplitude,

v_{v_{t}}

is the Gaussian white noise, and “

\circ

” denotes the Hadamard product.

As depicted in Figure 7, the simulation flight path primarily consisted of a straight flight and turns; the px4 autopilot controlled the SUAV to fly along the path at a cruise speed of 30 m/s under different wind field conditions. Data from the red route segments were intercepted for analysis.

4.3.2. Simulation Process

In real-world scenarios, SUAVs commonly employ autopilots to ensure safety and efficient flight operations in line with the mission. Consequently, the flow angles are often constrained within a narrow range, making it challenging to excite the dynamics of the SUAV. To replicate realistic conditions in this study, the “Mission” mode of the px4 autopilot was utilized for verification purposes. The attitude variations of the SUAV in the five scenarios are depicted in Figure 8.

During the simulation implementation, we also collected

[u, v, w]

and

[α, β]

from JSBSim, which jointly, with the preset aerodynamic derivatives, served as the reference for validating the OAFE method.

4.4. Simulation Results and Analysis

The OAFE method was evaluated by comparing the estimated values of

v_{r}^{b}

, the coefficient of side force

C_{Y}

, the coefficient of lift

C_{L}

{C_{L}}_{α}, C_{Y_{β}}

, and

f^{b, y}, f^{b, z}

(labeled as “OAFE”) with the reference values provided by the HIL system (labeled as “truth”). In the simulations, the OAFE output frequency was set to 20 Hz. Data from two complete turns along the above flight path were collected to assess the performance of the OAFE method. All the estimated values were highly accurate (Figure 9, Figure 10 and Figure 11).

As mentioned in Section 3.2, the air data parameters (including

v_{r}^{b}

) and aerodynamic parameters determine its aerodynamic force. A comparison was made between the calculated values of

f^{b, y}

and

f^{b, z}

(obtained from the estimated

[u, v, w]

and aerodynamic parameters) with the truth values. This comparison is illustrated in Figure 12.

As shown in Table 5, the OAFE method demonstrates a strong efficacy in estimating

v

w

β

, and

α

, across five distinct wind field conditions. The root mean square error (RMSE) of

α

varies less under different wind field conditions, whereas the accuracy of

β

is higher during flights without crosswinds.

C_{L_{α}}

and

C_{Y_{β}}

are theoretically constant values when

α

and

β

are small. Thus, the real-time estimates of

C_{L_{α}}

and

C_{Y_{β}}

obtained from the OAFE should converge to their true values, respectively. More importantly, in this study, the convergence of

C_{L_{α}}

and

C_{Y_{β}}

with different initial values can be a factor in demonstrating the feasibility of the OAFE. The simulation results (shown in Figure 13) demonstrate that

C_{L_{α}}

shows good convergence with different initial values in C5.

5. Field Tests

The ideas for evaluation in the flight tests are shown in Table 6. In the flight tests, the model-aided method provides reference values for

α

and

β

, as it is an extensively researched method with an accuracy within 0.5° and 1° [5,20,23]. What we aim to demonstrate in the field tests is that the proposed method achieves the estimation of flow angles without the prior acquisition of aerodynamic parameters, rather than achieving superior estimation accuracy.

Before conducting the field test, an aerodynamic derivatives model was developed for the test SUAV [37]. This model serves two purposes in the field tests, as follows: to provide reference aerodynamic parameters; and to complete a model-aided estimation. Appendix A.3 provides detailed aerodynamic parameters.

5.1. SUAV for Flight Tests

The Sula89 SUAV, employed for conducting the field tests, is presented in Figure 14. The Sula89 adopts a tandem layout with a high aspect ratio and equal-chord straight wings. The vertical tails act as a lateral stabilizer. Detailed geometric characteristics of the Sula89 are shown in Table 7.

5.2. Hardware Modules

As shown in Figure 14, the Sula89 is equipped with a Pixhawk2.1 autopilot (running px4 v1.8.0) that facilitates automatic flight control, and an NVIDIA Jetson Xavier NX3 that runs the OAFE. The decision not to integrate the OAFE with the EKF module in px4 was made to maintain the original navigation module, thereby ensuring flight safety, and minimizing dimensionality and computational burden on the estimator. The data used in this study were collected by the Pixhawk2.1, which was connected with a CUAV digital airspeed and pitot tube, as well as a HEX Here3, which provided satellite positioning information.

The IMU cube of the Pixhawk 2.1 was replaced with the 6-axis IMU (BMI160), as well as the STMicroelectronics 3-axis magnetometer LIS3MDL. The BMI160 is a 6-axis IMU consisting of a digital, tri-axial, 16-bit acceleration sensor and a tri-axial 16-bit gyroscope.

Table 3 in Section 4 illustrates the sensor data that can be obtained from the hardware modules. [

ϕ, θ, ψ

] were estimated based on the fusion of IMU, magnetometer and GPS data with high accuracy [38]. The sensor data and estimated [

ϕ, θ, ψ

] transmitted from px4 to NX3 were downsampled to 50 Hz. In the flight test conducted for this study, the Sula89 was not equipped with an air data system or a flow angle sensor.

5.3. Field Test Setup

In 2024, a field flight verification of the OAFE was conducted in Anhui Province, China, at a standard airplane flight field with legal airspace. The site is located on a huge farm, as indicated in Figure 15. The local temperature was 17 °C, and the weather was overcast with a level 2 east wind (about 3 m/s).

In the flight test, the initialization of the OAFE was consistent with the simulation, as in Section 4.1.

First, the structure and electrical components of the Sula89a were checked. Then, the SUAV took off remotely in stabilized mode, and subsequently switched to mission mode to fly autonomously along a preset four-sided route (shown in Figure 15). Finally, the SUAV was remotely controlled to land, after which the data stored onboard were archived. The specific control parameters of the px4 employed in this test can be found in [39].

5.4. Results and Analysis

The estimated

[u, v, w]

obtained from the OAFE method were compared with the reference, as depicted in Figure 16. The model-free-1 method, which is kinematic-based and assumes no wind (

v_{w}^{n} = 0

), is labeled with “mf1” in the figures below [16]. The model-free-2 method, with the assumption of

β = 0

and 2D horizontal wind, is labeled with “mf2” [15]. These two methods were selected based on the principle of immediate availability, like the OAFE. That is, there is no need for direct flow angle measurements, the prior acquisition of aerodynamic parameters, or pre-training.

The pitot airspeed is labeled “Pitot”. The

v_{r}^{b}

estimated based on the model-aided method is labeled as “model-aided” and was used as the reference. Consistent with the approach outlined in [25], we also utilized the pre-identified aerodynamic parameters of the Sula89 [39].

As shown in Figure 16, the mf1 method presents the largest error for the estimation of

[u, v, w]

. The error of

u

exhibits periodic oscillations due to the Sula89 flying a quadrilateral course in the wind field. The mf2 method shows significantly improved accuracy in estimating

u

, but not in estimating

v

and

w

. The OAFE method consistently shows good agreement with the reference values for the estimation of

[u, v, w]

throughout the entire flight. The

u

estimated by OAFE is a fusion of inertial and differential pressure information, which is smoother than the reference

u_{p}

. We think that it is a beneficial effect, as

u

is a long-term state and changes slowly.

The flow angle

[α, β]

, calculated from the estimated

[u, v, w]

, is shown in Figure 17.

Taking the results obtained from the “model-aided” method as a reference, the statistical performance of estimated

[v, w]

and

[α, β]

are summarized in Table 8.

Regarding the mf1 method, under the assumption of no wind, even with

v

w

errors as low as 1–2 m/s, the method resulted in a significant flow angle error for the test SUAV. The error of

β

even reached 10° at 65 s, presumably due to gusts of wind. This is to be expected since this method is designed for high-speed general aviation airplanes. As for the mf2 method, the estimated

v

and

β

fluctuated around zero, respectively, because it utilized

β

= 0 as a measurement. Although its assumptions reduced the RMSE of

v

and

β

, it did not lead to accurate estimates. For the estimation of

α

, the improvement in mf2 over mf1 was very small.

The OAFE method significantly improved the estimation accuracy of

β

and

α

. Throughout the process, the OAFE estimates fit the reference values better and captured more nuanced dynamical changes. The estimation results of both

β

and

α

are close to the accuracy limits of the reference (

1 °

and

0.5 °

). Even though the lateral aerodynamics adopt only a “

β - C_{Y}

” model, the OAFE still demonstrated good accuracy in estimating

β

. This is attributed to the OAFE’s ability to capture the rapid changes in

β

by fusing inertial data, even without incorporating aileron inputs, roll rates, or yaw rates into its aerodynamic model.

The distribution of estimation errors in the

α

domain is illustrated in Figure 18, revealing an increase in error with the increase in

α

. This may be attributed to the fact that the OAFE method utilizes a “

L - α

” linear model, while the tested SUAV’s “

L - α

” curve has some anticipated nonlinear characteristics.

The estimated

C_{Y}, C_{L}

by the OAFE are presented in Figure 19. The coefficients obtained from the CFD work were used as a reference and labeled as “CFD”. Among these, the lift coefficient

C_{L}

showed a more obvious hysteresis in the interval of

α = [0.05,0.12]

The estimated

f^{b, y}

and

f^{b, z}

, derived from the estimated

v_{r}^{b}

and estimated

[{C_{L}}_{0}, {C_{L}}_{α}, C_{L_{q}}, C_{L_{δ_{e}}}, C_{Y_{β}}]

, are depicted in Figure 20. The measurement components from the IMU serve as the reference for comparison.

For further analysis, the original sensor data stored in the TF memory card were utilized for offline validation. As shown in Figure 21, the results indicate that

[C_{Y_{β}}, C_{L_{α}}]

converge rapidly when different initial values are set. The amplitude of

C_{Y_{β}}

reduces to 7% after 3 s. Similarly, the amplitude of

C_{L_{α}}

decreases to 5% after 3 s.

6. Conclusions

An OAFE method that yields flow angle estimates of SUAVs solely using a standard sensor suite is proposed in this study. The OAFE’s filter loop utilizes inertial data and incorporates noisy aerodynamic derivatives identified online to refine its performance. The noisy aerodynamic derivatives are generated in real time by the optimizer loop of the OAFE in the absence of flow angle measurements. The optimizer approximately determines these derivatives based on multi-period sensor data, motion-continuity rules, and aerodynamics.

The OAFE method was evaluated through simulations and field tests. The simulation results revealed that the estimated

α

and

β

displayed an RMSE of approximately

0.11 °

and

0.24 °

, respectively. In field flight tests, the results indicated an RMSE of 0.75° for

β

and 0.59° for

α

, which fell within the error range of the reference method.

The OAFE method does not require flow angle measurements or the prior acquisition of aerodynamic parameters (derivatives or neural networks). Therefore, it is economical and does not require identification or training efforts in advance, making it suitable for quick deployment on different SUAVs and adaptable to changes in aerodynamic parameters caused by icing or other factors.

Further efforts will be dedicated to completing the estimation of other key aerodynamic derivatives without the need for α and β measurements. In addition, the OAFE will be utilized for the estimation of complete-state air data.

Author Contributions

Conceptualization, J.L. (Jie Li); Methodology, Z.W.; Software, Z.W. and Y.Y. (Yu Yang); Validation, X.W.; Formal analysis, J.L. (Juan Li); Investigation, C.L. and B.Y.; Resources, J.L. (Jie Li); Data curation, Y.Y. (Yu Yang) and X.W.; Writing—original draft, Z.W.; Writing—review & editing, J.L. (Juan Li) and Y.Y. (Yachao Yang); Supervision, C.L.; Project administration, C.L., Y.Y. (Yachao Yang) and B.Y.; Funding acquisition, J.L. (Juan Li). All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported in part by the National Natural Science Foundation of China under Grant 62373053.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Nomenclature

$v_{r}$	Relative velocity
$V_{a}$	Airspeed
$u_{p}$	Pitot airspeed
$U_{0}$	Trim airspeed
$α$	Angle of attack
$β$	Side slip angle
$\bar{q}$	Dynamic pressure
$v_{g}$	Ground velocity
$v_{w}$	Wind velocity over ground
$F_{a}$	Aerodynamic force vector
$ρ$	Air density
$\bar{c}$	Mean aerodynamic chord
$m$	Mass of SUAV
$\bar{S}$	Reference area
$G$	Gravity
$F_{T}$	Thrust
$C_{a}^{b}$	Rotation matrix from $a$ frame to $b$ frame
$T$	Sample time interval
$u$	Control input
$f (\cdot)$	Nonlinear system process equations
$h (\cdot)$	Measurement equations
$Q_{k - 1}$	Process error covariance matrix
$R_{k}$	Measurement error covariance matrix
$[δ_{e}, δ_{r}]$	Elevator and rudder deflections
$[L, D, Y]$	Lift, drag, and lateral forces
$[u, v, w]$	Relative velocity components along the x, y, and z axes in body frame
$ξ = {[ϕ, θ, ψ]}^{T}$	Euler angles
$g = {[g, 0, 0]}^{T}$	Gravitational acceleration
$ω_{i b}$ , $ω$ , or ${[p, q, r]}^{T}$	Angular rates
$f^{b} = {[f^{b, x}, f^{b, y}, f^{b, z}]}^{T}$	Specific force in body frame
$[C_{Y_{β}}, {C_{Y}}_{δ_{r}}, C_{L_{0}}, C_{L_{α}}, C_{L_{q}}, C_{L_{δ_{e}}}]$	Aircraft side and lift force coefficient derivatives
Superscripts:
$b$	denotes values in body frame
$n$	denotes values in navigation frame
$a$	denotes values in wind frame

Appendix A

Appendix A.1. Rotation Matrix

The conversion from wind frame to body frame is provided by the direction cosine matrix, which is defined by successive rotations of

α

and

β

as follows:

C_{a}^{b} = [\begin{matrix} \cos α & 0 & - \sin α \\ 0 & 1 & 0 \\ \sin α & 0 & \cos α \end{matrix}] [\begin{matrix} \cos β & - \sin β & 0 \\ \sin β & \cos β & 0 \\ 0 & 0 & 1 \end{matrix}]

(A1)

The conversion from navigation frame to body frame is provided by the direction cosine matrix, which is defined by successive rotations of the roll, pitch, and yaw angles

[φ, θ, ψ]

of the aircraft, as follows:

C_{n}^{b} = [\begin{matrix} 1 & 0 & 0 \\ 0 & \cos ϕ & \sin ϕ \\ 0 & - \sin ϕ & \cos ϕ \end{matrix}] [\begin{matrix} \cos θ & 0 & - \sin θ \\ 0 & 1 & 0 \\ \sin θ & 0 & \cos θ \end{matrix}] [\begin{matrix} \cos ψ & \sin ψ & 0 \\ - \sin ψ & \cos ψ & 0 \\ 0 & 0 & 1 \end{matrix}]

(A2)

Appendix A.2. Derivation of Process Noise of $v_{r}^{b}$

a. Derivation of Equation (35)

\begin{array}{l} Q_{v, k - 1}^{A} = E [\int_{t_{k - 1}}^{t_{k}} (δ_{ω^{b}} \times) d t \cdot {v_{r}^{b}}_{k - 1}] {[({v_{r}^{b}}_{k - 1}^{T}) \int_{t_{k - 1}}^{t_{k}} (δ_{ω^{b}} \times) d t]}^{T} \\ = E \{({v_{r}^{b}}_{k - 1} \times) \int_{t_{k - 1}}^{t_{k}} δ_{ω^{b}} d t {[\int_{t_{k - 1}}^{t_{k}} δ_{ω^{b}} d t]}^{T} {({v_{r}^{b}}_{k - 1} \times)}^{T}\} \\ = ({v_{r}^{b}}_{k - 1} \times) E \{\int_{t_{k - 1}}^{t_{k}} δ_{ω^{b}} d t {[\int_{t_{k - 1}}^{t_{k}} δ_{ω^{b}} d t]}^{T}\} {({v_{r}^{b}}_{k - 1} \times)}^{T} \end{array}

(A3)

Since

δ_{ω^{b}}

can be approximated as Gaussian white noise in the time domain and the 3-axis components are uncorrelated with each other, thus:

\{\begin{array}{l} E [δ_{ω^{b}} (t) \int_{τ_{k - 1}}^{τ_{k}} {δ_{ω^{b}}}^{T} (τ) d τ] = 0, i f : τ \neq t \\ E [δ_{ω^{b}} (t) \int_{τ_{k - 1}}^{τ_{k}} {δ_{ω^{b}}}^{T} (τ) d τ] = δ_{ω^{b}} (t) {δ_{ω^{b}}}^{T} (t), i f : τ \neq t \end{array}

(A4)

Therefore:

E \{\int_{t_{k - 1}}^{t_{k}} δ_{ω^{b}} d t {[\int_{t_{k - 1}}^{t_{k}} δ_{ω^{b}} d t]}^{T}\} = \int_{t_{k - 1}}^{t_{k}} δ_{ω^{b}} {δ_{ω^{b}}}^{T} d t

(A5)

Finally:

Q_{v, k - 1}^{A} = ({v_{r}^{b}}_{k - 1} \times) σ_{ω^{b}}^{2} {({v_{r}^{b}}_{k - 1} \times)}^{T} \cdot T

(A6)

b. Derivation of Equation (36)

\begin{array}{l} Q_{v, k - 1}^{C} = E \{\int_{t_{k - 1}}^{t_{k}} δ (C_{n}^{b} g^{n}) d t\} {\{\int_{t_{k - 1}}^{t_{k}} δ (C_{n}^{b} g^{n}) d t\}}^{T} \\ = E \{\int_{t_{k - 1}}^{t_{k}} C_{n}^{b} (δ_{ξ} \times) g^{n} d t\} {\{\int_{t_{k - 1}}^{t_{k}} C_{n}^{b} (δ_{ξ} \times) g^{n} d t\}}^{T} \\ = E \{C_{n}^{b} (g^{n} \times) \int_{t_{k - 1}}^{t_{k}} δ_{ξ} d t\} \{\int_{t_{k - 1}}^{t_{k}} {δ_{ξ}}^{T} d t (g^{n} \times) C_{b}^{n}\} \end{array}

(A7)

where

δ_{ξ}

denotes the attitude error.

The same principle as a is applicable, as follows:

Q_{v, k - 1}^{C} = C_{n}^{b} (g^{n} \times) σ_{ξ}^{2} (g^{n} \times) C_{b}^{n} \cdot T

(A8)

Appendix A.3. Aerodynamic Parameters

Table A1. Aerodynamic force and moment derivatives.

Components	Values	Components	Values	Components	Values	Components	Values
$C_{L_{0}}$	0.52	$C_{Y_{p}}$	−0.20	$C_{l_{δ_{a}}}$	0.19	$C_{m_{δ_{e}}}$	−4.44
${C_{L}}_{α}$	5.10	$C_{D_{0}}$	0.07	$C_{m_{0}}$	−0.23	$C_{m_{\dot{α}}}$	−6.97
$C_{L_{q}}$	4.67	$C_{D_{α^{2}}}$	3.31	$C_{m_{α}}$	−0.90	$C_{n_{β}}$	0.37
$C_{L_{δ_{e}}}$	1.07	$C_{D_{β}}$	0.17	$C_{m_{α^{2}}}$	−6.50	$C_{n_{r}}$	−0.003
$C_{Y_{β}}$	−0.53	$C_{l_{β}}$	−0.14	$C_{m_{α^{3}}}$	−35.89	$C_{n_{δ_{r}}}$	−0.06
$C_{Y_{r}}$	0.14	$C_{l_{p}}$	−0.59	$C_{m_{q}}$	−89.23

α \in [- 0.14, 0.21] .

References

Kai, J.-M.; Hamel, T.; Samson, C. A Unified Approach to Fixed-Wing Aircraft Path Following Guidance and Control. Automatica 2019, 108, 108491. [Google Scholar] [CrossRef]
Yang, D.; Zang, J.; Liu, J.; Liu, K. Time-Domain Identification Method Based on Data-Driven Intelligent Correction of Aerodynamic Parameters of Fixed-Wing UAV. Drones 2023, 7, 594. [Google Scholar] [CrossRef]
Liu, F.; Lu, L.; Zhang, Z.; Xie, Y.; Chen, J. Intelligent Trajectory Prediction Algorithm for Hypersonic Vehicle Based on Sparse Associative Structure Model. Drones 2024, 8, 505. [Google Scholar] [CrossRef]
Sankaralingam, L.; Ramprasadh, C. A Comprehensive Survey on the Methods of Angle of Attack Measurement and Estimation in UAVs. Chin. J. Aeronaut. 2020, 33, 749–770. [Google Scholar] [CrossRef]
Tian, P.; Chao, H.; Flanagan, H.P.; Hagerott, S.G.; Gu, Y. Design and Evaluation of UAV Flow Angle Estimation Filters. IEEE Trans. Aerosp. Electron. Syst. 2019, 55, 371–383. [Google Scholar] [CrossRef]
Sankaralingam, L.; Ramprasadh, C. Angle of Attack Measurement Using Low-Cost 3D Printed Five Hole Probe for UAV Applications. Measurement 2021, 168, 108379. [Google Scholar] [CrossRef]
Borup, K.T.; Fossen, T.I.; Johansen, T.A. A Machine Learning Approach for Estimating Air Data Parameters of Small Fixed-Wing UAVs Using Distributed Pressure Sensors. IEEE Trans. Aerosp. Electron. Syst. 2020, 56, 2157–2173. [Google Scholar] [CrossRef]
Mersha, B.W.; Jansen, D.N.; Ma, H. Angle of Attack Prediction Using Recurrent Neural Networks in Flight Conditions with Faulty Sensors in the Case of F-16 Fighter Jet. Complex Intell. Syst. 2023, 9, 2599–2611. [Google Scholar] [CrossRef]
Liu, Y.; Zhang, C.; Yan, X.; Liu, W. Flush Air Data Sensing Based on Dimensionless Input and Output Neural Networks with Less Data. IEEE Trans. Aerosp. Electron. Syst. 2022, 59, 1411–1425. [Google Scholar] [CrossRef]
Popowski, S.; Dąbrowski, W. Measurement and Estimation of the Angle of Attack and the Angle of Sideslip. Aviation 2015, 19, 19–24. [Google Scholar] [CrossRef]
Lim, H.; Ryu, H.; Rhudy, M.B.; Lee, D.; Jang, D.; Lee, C.; Park, Y.; Youn, W.; Myung, H. Deep Learning-Aided Synthetic Airspeed Estimation of UAVs for Analytical Redundancy With a Temporal Convolutional Network. IEEE Robot. Autom. Lett. 2022, 7, 17–24. [Google Scholar] [CrossRef]
Wenz, A.; Johansen, T.A.; Cristofaro, A. Combining Model-Free and Model-Based Angle of Attack Estimation for Small Fixed-Wing UAVs Using a Standard Sensor Suite. In Proceedings of the 2016 International Conference on Unmanned Aircraft Systems (ICUAS), Arlington, VA, USA, 7–10 June 2016; pp. 624–632. [Google Scholar]
Wenz, A.; Johansen, T.A. Moving Horizon Estimation of Air Data Parameters for UAVs. IEEE Trans. Aerosp. Electron. Syst. 2020, 56, 2101–2121. [Google Scholar] [CrossRef]
Cho, A.; Kim, J.; Lee, S.; Kee, C. Wind Estimation and Airspeed Calibration Using a UAV with a Single-Antenna GPS Receiver and Pitot Tube. IEEE Trans. Aerosp. Electron. Syst. 2011, 47, 109–117. [Google Scholar] [CrossRef]
Yang, Y.; Liu, X.; Liu, X.; Guo, Y.; Zhang, W. Model-Free Integrated Navigation of Small Fixed-Wing UAVs Full State Estimation in Wind Disturbance. IEEE Sens. J. 2022, 22, 2771–2781. [Google Scholar] [CrossRef]
Lerro, A.; Brandl, A.; Gili, P. Model-Free Scheme for Angle-of-Attack and Angle-of-Sideslip Estimation. J. Guid. Control Dyn. 2021, 44, 595–600. [Google Scholar] [CrossRef]
Bagherzadeh, S.A. Nonlinear Aircraft System Identification Using Artificial Neural Networks Enhanced by Empirical Mode Decomposition. Aerosp. Sci. Technol. 2018, 75, 155–171. [Google Scholar] [CrossRef]
Freeman, D.B. Angle of Attack Computation System; AFFDL-TR-73-89; Air Force Flight Dynamics Laboratory: Wright-Patterson AFB, OH, USA, 1973. [Google Scholar]
Tian, P.; Chao, H. Model Aided Estimation of Angle of Attack, Sideslip Angle, and 3D Wind without Flow Angle Measurements. In Proceedings of the 2018 AIAA Guidance, Navigation, and Control Conference, Kissimmee, FL, USA, 8–12 January 2018. [Google Scholar]
Youn, W.; Choi, H.; Cho, A.; Kim, S.; Rhudy, M.B. Aerodynamic Model-Aided Estimation of Attitude, 3-D Wind, Airspeed, AOA, and SSA for High-Altitude Long-Endurance UAV. IEEE Trans. Aerosp. Electron. Syst. 2020, 56, 4300–4314. [Google Scholar] [CrossRef]
Youn, W.; Lim, H.; Choi, H.S.; Rhudy, M.B.; Ryu, H.; Kim, S.; Myung, H. State Estimation for HALE UAVs With Deep-Learning-Aided Virtual AOA/SSA Sensors for Analytical Redundancy. IEEE Robot. Autom. Lett. 2021, 6, 5276–5283. [Google Scholar] [CrossRef]
Karali, H.; Uzun, M.; Yuksek, B.; Inalhan, G. Data-Driven Synthetic Air Data Estimation System Development for a Fighter Aircraft. In Proceedings of the AIAA AVIATION 2023 Forum, American Institute of Aeronautics and Astronautics. San Diego, CA, USA, 12–16 June 2023; p. 3439. [Google Scholar]
Youn, W.; Choi, H.S.; Ryu, H.; Kim, S.; Rhudy, M.B. Model-Aided State Estimation of HALE UAV With Synthetic AOA/SSA for Analytical Redundancy. IEEE Sens. J. 2020, 20, 7929–7940. [Google Scholar] [CrossRef]
Lu, H.; Gao, L.; Yan, Y.; Hou, M.; Wang, C. Wind Disturbance Compensated Path-Following Control for Fixed-Wing UAVs in Arbitrarily Strong Winds. Chin. J. Aeronaut. 2023, 37, 431–445. [Google Scholar] [CrossRef]
Langelaan, J.W.; Alley, N.; Neidhoefer, J. Wind Field Estimation for Small Unmanned Aerial Vehicles. J. Guid. Control Dyn. 2011, 34, 1016–1030. [Google Scholar] [CrossRef]
PX4 Autopilot Software. Available online: https://github.com/PX4/PX4-Autopilot (accessed on 1 July 2023).
Yuan, L.; Zheng, J.; Wang, X.; Ma, L. Attitude Control of a Mass-Actuated Fixed-Wing UAV Based on Adaptive Global Fast Terminal Sliding Mode Control. Drones 2024, 8, 305. [Google Scholar] [CrossRef]
Seo, G.; Kim, Y.; Saderla, S. Kalman-Filter Based Online System Identification of Fixed-Wing Aircraft in Upset Condition. Aerosp. Sci. Technol. 2019, 89, 307–317. [Google Scholar] [CrossRef]
Arasaratnam, I.; Haykin, S. Cubature Kalman Filters. IEEE Trans. Autom. Control 2009, 54, 1254–1269. [Google Scholar] [CrossRef]
Li, K.; Chen, X.; Liu, H.; Wang, S.; Li, K.; Li, B. Performance Analysis of the Thermal Automatic Tracking Method Based on the Model of the UAV Dynamic Model in a Thermal and Cubature Kalman Filter. Drones 2023, 7, 102. [Google Scholar] [CrossRef]
Beard, R.W.; McLain, T.W. Small Unmanned Aircraft; Princeton University Press: Princeton, NJ, USA, 2012; ISBN 978-0-691-14921-9. [Google Scholar]
Bonnor, N. Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems, 2nd ed.; Groves, P.D., Ed.; Artech House: London, UK, 2013; p. 776. ISBN 978-1-60807-005-3. [Google Scholar]
Wenz, A.; Johansen, T.A. Estimation of Wind Velocities and Aerodynamic Coefficients for UAVs Using Standard Autopilot Sensors and a Moving Horizon Estimator. In Proceedings of the 2017 International Conference on Unmanned Aircraft Systems (ICUAS), Miami, FL, USA, 13–16 June 2017; pp. 1267–1276. [Google Scholar]
Minaev, E.; Quijada Pioquinto, J.G.; Shakhov, V.; Kurkin, E.; Lukyanov, O. Airfoil Optimization Using Deep Learning Models and Evolutionary Algorithms for the Case Large-Endurance UAVs Design. Drones 2024, 8, 570. [Google Scholar] [CrossRef]
Forster, C.; Carlone, L.; Dellaert, F.; Scaramuzza, D. On-Manifold Preintegration for Real-Time Visual--Inertial Odometry. IEEE Trans. Robot. 2017, 33, 1–21. [Google Scholar] [CrossRef]
Nijboer, J.; Armanini, S.F.; Karasek, M.; de Visser, C.C. Longitudinal Grey-Box Model Identification of a Tailless Flapping-Wing MAV Based on Free-Flight Data. In Proceedings of the AIAA Scitech 2020 Forum; American Institute of Aeronautics and Astronautics, Orlando, FL, USA, 6–10 January 2020. [Google Scholar]
Wang, Z.; Xiong, J.; Cheng, X.; Li, J. Establishment and Verification of Longitudinal Aerodynamic Model of Tandem Wing Aircraft. IOP Conf. Ser. Mater. Sci. Eng. 2019, 563, 032022. [Google Scholar] [CrossRef]
Zhao, Y. Performance Evaluation of Cubature Kalman Filter in a GPS/IMU Tightly-Coupled Navigation System. Signal Process. 2016, 119, 67–79. [Google Scholar] [CrossRef]
Yang, Y.; Liu, C.; Li, J.; Yang, Y.; Li, J.; Zhang, Z.; Ye, B. Design, Implementation, and Verification of a Low-cost Terminal Guidance System for Small Fixed-wing UAVs. J. Field Robot. 2021, 38, 801–827. [Google Scholar] [CrossRef]

Figure 1. Illustration of coordinates and flow angles (

D

Y

, and

L

are the drag, lateral, and lift forces).

Figure 1. Illustration of coordinates and flow angles (

D

Y

, and

L

are the drag, lateral, and lift forces).

Figure 2. OAFE method structure.

Figure 3. Pseudo-measurements of [

{\tilde{C}}_{Y_{β}} {\tilde{C}}_{L}_{α}

] in flight tests.

Figure 3. Pseudo-measurements of [

{\tilde{C}}_{Y_{β}} {\tilde{C}}_{L}_{α}

] in flight tests.

Figure 4. The construction of the unknown sequence

{\tilde{v}}_{r, i}^{b} .

Yellow arrows are unused recursive constructors, while the blue arrows are the direct constructors.

Figure 4. The construction of the unknown sequence

{\tilde{v}}_{r, i}^{b} .

Yellow arrows are unused recursive constructors, while the blue arrows are the direct constructors.

Figure 5. Hardware scheme for the simulation.

Figure 6. Simulation framework.

Figure 7. Waypoints for simulation setup.

Figure 8.

ϕ

and

θ

in different wind fields.

Figure 8.

ϕ

and

θ

in different wind fields.

Figure 9. Comparison of estimated

v

and

w

with the reference.

Figure 9. Comparison of estimated

v

and

w

with the reference.

Figure 10. Comparison of estimated

β

and

α

with the reference.

Figure 10. Comparison of estimated

β

and

α

with the reference.

Figure 11. Comparison of estimated

C_{Y}

and

C_{L}

with the reference.

Figure 11. Comparison of estimated

C_{Y}

and

C_{L}

with the reference.

Figure 12. Comparison of estimated

f^{b, y}

and

f^{b, z}

with the reference.

Figure 12. Comparison of estimated

f^{b, y}

and

f^{b, z}

with the reference.

Figure 13. Convergence of

C_{L_{α}}

with different initial values.

Figure 13. Convergence of

C_{L_{α}}

with different initial values.

Figure 14. Fixed-wing SUAV for field tests.

Figure 15. (a) standard model airplane flight field; and (b) four-side route in flight test.

Figure 16. Comparison of

[u, v, w]

results from OAFE with the reference in the flight test.

Figure 16. Comparison of

[u, v, w]

results from OAFE with the reference in the flight test.

Figure 17. Comparison of

[α, β]

results from the OAFE with the reference in flight tests.

Figure 17. Comparison of

[α, β]

results from the OAFE with the reference in flight tests.

Figure 18. “Error-

α

” relationship of estimated

w

. The orange straight line is a linear fit between the error of

w

and

α

Figure 18. “Error-

α

” relationship of estimated

w

. The orange straight line is a linear fit between the error of

w

and

α

Figure 19. Estimated

C_{Y}

and

C_{L}

compared with CFD results.

Figure 19. Estimated

C_{Y}

and

C_{L}

compared with CFD results.

Figure 20. Estimated

f^{b, y}

and

f^{b, z}

compared with measurements.

Figure 20. Estimated

f^{b, y}

and

f^{b, z}

compared with measurements.

Figure 21. Convergence of

C_{Y_{β}}

and

C_{L_{α}}

with different initial values.

Figure 21. Convergence of

C_{Y_{β}}

and

C_{L_{α}}

with different initial values.

Table 1. Validation methods.

OAFE Estimates	Reference for Simulation
$u, v, w, α, β, f^{b}$	Truth values from the HIL Syst
$C_{L}, C_{β}, {C_{L}}_{α}, C_{Y_{β}}$	Derivatives of SUAV in the HIL Syst

Table 2. Noise characteristics of the OAFE.

Input Noise $σ^{2}$		$Measurement Noise σ^{2}$
$p$	${0.00109}^{2} {(r a d / s)}^{2}$	$u_{p}$	${0.2}^{2} {(m / s)}^{2}$
$q$	${0.00109}^{2} {(r a d / s)}^{2}$	$f^{b, y}$	$4 \times 10^{- 5} {(m / s^{2})}^{2}$
$r$	${0.00109}^{2} {(r a d / s)}^{2}$	$f^{b, z}$	$4 \times 10^{- 5} {(m / s^{2})}^{2}$
$ϕ$	${0.014}^{2} {r a d}^{2}$	${\tilde{C}}_{L}_{α}$	$R_{p, 11}$ (Equation (50))
$θ$	${0.014}^{2} {r a d}^{2}$	${\tilde{C}}_{Y_{β}}$	$R_{p, 22}$ (Equation (50))
$ψ$	${0.009}^{2} {r a d}^{2}$
$δ_{e}$	$8 \times 10^{- 5} r a d^{2}$
$δ_{r}$	$8 \times 10^{- 5} r a d^{2}$

Table 3. Sensor’s performance parameters.

Sensors	Type	Noise/Error
BMI160	IMU	$0.014 ° / S / \sqrt{H z}$
BMI160	IMU	$150 μ g / \sqrt{H z}$
MS5525	pressure sensor and pitot tube	$\pm 0.84 P a$
HEX Here3	GPS Module	$pos : 2.5 m$ $vel : 0.1 m / s$

Table 4. Wind fields set-up in flight simulation.

Wind Field Conditions	Wind Fields (m/s)
Case 1	$[0, 0, 0]$
Case 2	$[0, 3, 0]$
Case 3	$[0, 0, - 3]$
Case 4	$[0, 2.2, - 2.2]$
Case 5	$[- 6, 0, 0]$

Table 5. Root mean square errors of estimates.

Case	$v$ (m/s)	$w$ (m/s)	$β$ (°)	$α$ (°)
C1	0.038	0.127	0.072	0.239
C2	0.041	0.151	0.076	0.234
C3	0.028	0.137	0.051	0.257
C4	0.076	0.091	0.141	0.170
C5	0.069	0.155	0.129	0.293

Table 6. Validation methods for flight tests.

OAFE Estimates	Reference for Field Test
$u, v, w$ $α, β$	Model-aided method
$C_{L}, C_{β}, {C_{L}}_{α}, C_{Y_{β}}$	Parameters from a CFD work; Derivatives convergence
$f^{b, y}, f^{b, z}$	$f^{b, y}, f^{b, z}$ measured by IMU

Table 7. Geometric characteristics of the Sula89.

Parameters	Values
Aerodynamic configuration	Tandem wing
Fuselage length	0.92 m
Wingspan	1.56 m
Reference area	$0.26 m^{2}$
Aspect ratio	11
Mean aerodynamic chord	0.19 m
Mass	8.9 kg
Cruise speed	30 m/s

Table 8. RMSE of estimates from different methods.

Estimates	$v$ (m/s)	$w$ (m/s)	$β$ (°)	$α$ (°)
mf1	1.920	1.062	3.581	1.909
mf2	0.987	1.001	1.830	1.837
OAFE	0.407	0.324	0.748	0.593

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Wang, Z.; Li, J.; Liu, C.; Yang, Y.; Li, J.; Wu, X.; Yang, Y.; Ye, B. Optimization-Assisted Filter for Flow Angle Estimation of SUAV Without Adequate Measurement. Drones 2024, 8, 758. https://doi.org/10.3390/drones8120758

AMA Style

Wang Z, Li J, Liu C, Yang Y, Li J, Wu X, Yang Y, Ye B. Optimization-Assisted Filter for Flow Angle Estimation of SUAV Without Adequate Measurement. Drones. 2024; 8(12):758. https://doi.org/10.3390/drones8120758

Chicago/Turabian Style

Wang, Ziyi, Jie Li, Chang Liu, Yu Yang, Juan Li, Xueyong Wu, Yachao Yang, and Bobo Ye. 2024. "Optimization-Assisted Filter for Flow Angle Estimation of SUAV Without Adequate Measurement" Drones 8, no. 12: 758. https://doi.org/10.3390/drones8120758

APA Style

Wang, Z., Li, J., Liu, C., Yang, Y., Li, J., Wu, X., Yang, Y., & Ye, B. (2024). Optimization-Assisted Filter for Flow Angle Estimation of SUAV Without Adequate Measurement. Drones, 8(12), 758. https://doi.org/10.3390/drones8120758

Article Menu

Optimization-Assisted Filter for Flow Angle Estimation of SUAV Without Adequate Measurement

Abstract

1. Introduction

2. Problem Formulation

2.1. Kinematics Model

2.2. Dynamics Model

2.3. Dynamics Simplification

2.4. Cubature Kalman Filter (CKF)

3. OAFE Method

3.1. Filter Loop

3.1.1. State Dynamic and Measurement

3.1.2. Filtering Noise

3.2. Optimizer Loop

3.2.1. Objective Function

3.2.2. Improving Optimization Efficiency

3.2.3. Noises of Selected Derivatives

4. Simulation

4.1. Initialization of OAFE

4.2. HIL Simulation System

4.3. Simulation Setup

4.3.1. Simulation Scenes

4.3.2. Simulation Process

4.4. Simulation Results and Analysis

5. Field Tests

5.1. SUAV for Flight Tests

5.2. Hardware Modules

5.3. Field Test Setup

5.4. Results and Analysis

6. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

Nomenclature

Appendix A

Appendix A.1. Rotation Matrix

Appendix A.2. Derivation of Process Noise of v r b

Appendix A.3. Aerodynamic Parameters

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI

Appendix A.2. Derivation of Process Noise of $v_{r}^{b}$