Nomenclature

A

accelerometers output

a _x , a _y , a _z

acceleration components along body axes (m/sec²)

AD

air data

AOA

angle-of-attack (rad)

AOSS

Angle-of-sideslip (rad)

A _x , A _y , A _z

specific force components

EW-AEKF

exponentially weighted adaptive extended Kalman filter

f, h

non-linear functions

FDD

fault detection and diagnosis

FDI

fault detection and isolation

g

gravitational acceleration (m/sec²)

GPS

global positioning system

i

air data sensor index

IMU

inertial measurement unit

IRU

inertial reference unit

INS

inertial navigation system

J _T

threshold

k

time instant

p, q, r

body rotational rate components (rad/sec)

Q

process noise covariance matrix

R

measurement noise covariance matrix

PADS

primary air data sensors

RMSE

root mean square error

t

time (sec)

u _b, v _b, w _b

true airspeed velocity components (m/sec)

u

control input vector

U _g

ground speed (m/sec)

$U_{gx}^b$ , $U_{gy}^b$ , $\;U_{gz}^b$

ground speed velocity components along body axis (m/sec)

w & v

input and measurement noise vector

V _wx, V _wy, V _wz

wind velocity components in earth frame (m/sec)

U _t

true airspeed (m/sec)

V _w

wind speed (m/sec)

x & z

state and measurement vector

$\omega$

rate gyros output

$\alpha$

angle-of-attack sensor

$\beta$

angle-of-sideslip sensor

$\phi$ , $\theta$ , $\psi$ ;

attitude angles (rad)

e _Ut, $e_{\alpha}$ $e_{\beta}$

residual signals of air data sensors

$\mu$ & $\sigma$

mean and standard deviation of error

2.0 Problem statement

In general, model-based FDD techniques of ADSS are evaluated through the residual signals, which represent the inconsistency between the actual system behavior and the mathematical model of the system. For the sensor FDD of the proposed method, it can be defined as the difference between true air data sensor measurements and the estimated air data. The key advantages of using residuals are their ease of implementation and the fact that they carry exact information about faults such as stuck, bias, drift, increased noise faults, etc. It can be determined as,

(1) \begin{align}e(k) = {z_k} - {\hat z_k}\end{align}

where e (k) is the residual value at the k ^th time step, ${z_k}$ is the true sensor measurements, and ${\hat z_k}$ is the estimated measurement of the filter, such as an extended Kalman filter [Reference Julier and Uhlmann20]. The kinematic model-based primary air data sensors FDD (PADS-FDD) proposed in this paper consider the residuals of each sensor category of the air data system. Since the above Equation (1) is the generic form of residual signals, a complete knowledge of the signals and underlying dynamics is required to provide good cross-platform portability [Reference Guo, Zhong and Zhou21]. The residuals used in this FDD algorithm are presented as follows,

(2) \begin{align}\left\{ \begin{array}{l}{e_{{U_t}}} = {U_t} - {{\hat U}_t}\\ \\[-9pt]{e_\alpha } = \alpha - \hat \alpha \\ \\[-9pt]{e_\beta } = \beta - \hat \beta \end{array} \right.\end{align}

where U _t , $\alpha$ , $\beta$ are air data system components (true airspeed, AOA, and AOSS) and ${e_{{U_t}}}$ , ${e_\alpha }$ , ${e_\beta }$ are the residuals of true airspeed, AOA, and AOSS, respectively. The calculated residuals are used to monitor each sensor’s health status robustly. The air data is adaptively estimated, and residuals are generated by the adaptive version of the EKF, namely the EW-AEKF. The predetermined thresholds are fixed to detect various faults on the sensors. A fault occurs when residuals exceed the thresholds and are detected using the detection logic (discussed in Section 5).

In this study, we tested various fault scenarios for air data sensors that include freezing (stuck) fault, bias, drift and increased noise faults. To detect every fault efficiently, the threshold values are carefully determined, considering the sensor’s impact on the aircraft’s performance and control. In practice, the allowable error limit for sensors on civil aircraft is defined by civil aviation authorities such as the Civil Aviation Authority (CAA, UK), the European Aviation Safety Agency (EASA), the Federal Aviation Administration (FAA), etc. For example, the airspeed error limit that is allowable as per EASA’s Certification Specifications (CS 25.1323) of subpart F is ± 3% or ± 2.5m/s [22,23]. With this knowledge, the details of fixed thresholds are presented in Section 5 of this paper. In order to implement robust and reliable fault detection (PADS-FDD), we used the available sensor measurements other than the redundant air data (AD) sensor measurements. More clearly, the proposed AD estimation and FDD algorithm use measurements from the IRU and an integrated INS/GPS navigation system, which is less sensitive to external disturbances and readily available on the air data and inertial reference system (ADIRS) of modern aircraft. The following relationships are used to compute the air data from the estimated parameters. Using the velocity triangle relationship, one can extract the true airspeed as presented by,

(3) \begin{align}{U_g} = {U_t} + {V_w}\end{align}

where U _g is the ground speed measured in the earth frame, U _t is the true airspeed measured in the aircraft body frame, and V _w is the wind speed measured in the earth frame. The true airspeed can be written in component form by using airspeed and airflow angles as,

(4) \begin{align}\left\{ \begin{array}{l}u = {U_t}\cos \alpha \cos \beta \\ \\[-9pt]v = {U_t}\sin \beta \\ \\[-9pt]w = {U_t}\sin \alpha \cos \beta \end{array} \right.\end{align}

To extract the true airspeed without the aid of redundant Pitot tube measurements of the air data system, the Equation (3) is rewritten with a direction cosine matrix [Reference Stevens, Lewis and Johnson24] as,

(5) \begin{align}{U_t} = {\bf{C}}_e^b{U_g} - V_w^b\end{align}

where ${\bf{C}}_n^b$ is the transformation matrix from the earth frame to the body frame, and $V_w^b$ is the wind speed represented in the body frame. The airflow angles of the aircraft air data system are computed using the following relationship using body-axis velocity components:

(6) \begin{align}\left\{ \begin{array}{l}\alpha = {\tan ^{ - 1}}\left( {\dfrac{{{w_b}}}{{{u_b}}}} \right)\\ \\[-9pt]\beta = {\sin ^{ - 1}}\left( {\dfrac{{{v_b}}}{{{U_t}}}} \right)\end{array} \right.\end{align}

where $\alpha$ is the angle-of-attack, $\beta$ is the angle-of-side slip, and u _b , v _b , and w _b are the airspeed velocity components represented in the aircraft body frame. In Section 3, we will present the mathematical model for air data estimation.

3.0 Kinematic model and AD estimation

Kinematic models are indeed a special form of the general equations of motion (EOM) of an aircraft, for which specific forces (the outputs of accelerometers in the center of gravity) and angular rates serve as inputs. They can conveniently be written in state-space form. An aircraft EOM takes the form of three sets of first-order differential equations, one each for translational velocities, angular velocities, and attitude angles [Reference Mulder, Chu, Sridhar, Breeman and Labanm25]. For the air data estimation of the proposed work, we considered only the differential equations of velocities and attitude angles. The main advantage of a kinematic model is that it does not take inputs from flight control surface deflections (rudder, elevator and aileron) and is free from modeling uncertainties that degrade the model-based sensor FDD algorithms. Using the customary body-fixed reference frame, the equations for the components u, v, and w of true airspeed (U _t ) along the body axes take the following form:

(7) \begin{align}\left\{ \begin{array}{l}X = m(\dot u + qw - rv) + mg\sin \theta \\ \\[-8pt]Y = m(\dot v + ru - pw) - mg\cos \theta \sin \phi \\ \\[-8pt]Z = m(\dot w + pv - qu) - mg\cos \theta \cos \phi \end{array} \right.\end{align}

where p, q, and r denote the rotational rates about the aircraft body frame; $\theta$ and $\phi$ denote the pitch and roll angles, respectively; m denotes the aircraft mass; and g denotes the local acceleration due to gravity. X, Y and Z represent the components of the total aerodynamic forces, including the aerodynamic effects of propulsion systems. The aircraft kinematic-based model is driven by specific forces measured by ideal accelerometers and rotational rates measured by rate gyroscopes. In flight, such ideal accelerometers would measure the specific aerodynamic forces according to:

(8) \begin{align}\left\{ \begin{array}{l}X = {A_x}m\\ \\[-8pt]Y = {A_y}m\\ \\[-8pt]Z = {A_z}m\end{array} \right.\end{align}

where A _x , A _y , and A _z denote the specific aerodynamic forces along the aircraft body axes. On substituting Equation (8) into Equation (7) by eliminating the mass ‘m’ of the aircraft, we get

(9) \begin{align}\left\{ \begin{array}{l}\dot u = {a_x} + vr - wq - g\sin \theta \\ \\[-8pt]\dot v = {a_y} - ur + wp + g\cos \theta \sin \phi \\ \\[-8pt]\dot w = {a_z} + uq - vp + g\cos \theta \cos \phi \end{array} \right.\end{align}

The above set of equations represents the kinematical relations of aircraft velocities. The other set of equations is that governed by the first-order differential equations for the attitude (Euler) angles in the earth-fixed reference frame and is presented by

(10) \begin{align}\left\{ \begin{array}{l}\dot \phi = p + q\sin \phi \tan \theta + r\cos \phi \tan \theta \\ \\[-8pt]\dot \theta = q\cos \phi - r\sin \phi \\ \\[-8pt]\dot \psi = q\sin \phi \sec \theta + r\cos \phi \sec \theta \end{array} \right.\end{align}

The above Equations (9) and (10) can be called the kinematic model required for the estimation of air data. The development of a robust air data model is further enhanced with available sensor measurements, such as IRU and INS/GPS, to guarantee the reliability of estimation. Using Equation (9), the velocity dynamics required for the proposed method is written in kinematical form as

(11) \begin{align}\left\{ \begin{array}{l}\dot U_{gx}^b = {a_x} + U_{gy}^br - U_{gz}^bq - g\sin \theta \\ \\[-8pt]\dot U_{gy}^b = {a_y} - U_{gx}^br + U_{gz}^bp + g\cos \theta \sin \phi \\ \\[-8pt]\dot U_{gz}^b = {a_z} + U_{gx}^bq - U_{gy}^bp + g\cos \theta \cos \phi \end{array} \right.\end{align}

where $U_{gx}^b$ , $U_{gy}^b$ , $U_{gz}^b$ are the ground speed velocity components transformed from the earth frame to aircraft’s body frame which are readily available from the flight simulator sensor data set (in the earth frame). The above velocity dynamics include the effects of wind variations represented in the aircraft body frame. With the knowledge of Equations (10) and (11), the rest of the paper will focus on estimation of air data using the above kinematic relations. Note that the ground speed components U _gx , U _gy , U _gz are readily available from the flight sensor data of aircraft, either from flight simulation or real flight. Hence the output equations are presented as,

(12) \begin{align}\begin{array}{l}{U_t} = \sqrt {{U_{tx}}^2 + {U_{ty}}^2 + {U_{tz}}^2} \\ \\[-8pt]\alpha = {\tan ^{ - 1}}\left( {\dfrac{{{U_{tz}}}}{{{U_{tx}}}}} \right)\\ \\[-8pt]\beta = {\sin ^{ - 1}}\left( {\dfrac{{{U_{ty}}}}{{{U_t}}}} \right)\end{array}\end{align}

where U _tx , U _ty , and U _tz are the true airspeed components computed using the ground velocities and wind velocities in the body frame. Using Equations (4) and (11), the dynamics of air data variables [U _t , $\alpha$ , $\beta$ ] ^T i.e., $\dot U$ _t , $\dot \alpha $ , and $\dot \beta $ can be obtained as,

(13) \begin{align}{\dot U_t} = \dot U_{gx}^b\cos \alpha \cos \beta + \dot U_{gy}^b\sin \beta + \dot U_{gz}^b\sin \alpha \cos \beta \end{align}

(14) \begin{align}\dot \alpha = \frac{{\dot U_{gz}^b\cos \alpha - \dot U_{gx}^b\sin \alpha }}{{{U_t}\cos \beta }}\end{align}

(15) \begin{align}\dot \beta = \frac{{ - \dot U_{gx}^b\cos \alpha \sin \beta + \dot U_{gy}^b\cos \beta - \dot U_{gz}^b\sin \alpha \sin \beta }}{{{U_t}}}\end{align}

From the above information, the non-linear state-space model for the proposed AD estimation with the kinematics of the aircraft can be presented as,

(16) \begin{align}{\dot{\bf{x}}}\left( t \right) &= f\left[ {{\bf{x}}\left( t \right) + G\left( {{\bf{x}}\left( t \right)} \right)\,\left( {{\bf{u}}\left( t \right) + {\bf{w}}\left( t \right)} \right)} \right]\nonumber\\[4pt]{\bf{z}}\left( t \right) &= h\left( {{\bf{x}}\left( t \right)} \right)\nonumber\\[4pt]{{\bf{z}}_m}\left( t \right) &= \left[ {h\left( {{\bf{x}}\left( t \right)} \right) + {\bf{v}}\left( t \right)} \right]\end{align}

where x represents the state vector, u represents the input vector, z represents the measurement vector, and w and v represent the process and measurement noise vectors with zero-mean Gaussian with covariance matrices Q and R, respectively. f and h are the non-linear functions. The variables of system equations are written as follows:

(17) \begin{align}{\bf{x}} = {\left[ {{U_t},\alpha ,{\rm{ }}\beta ,{\rm{ }}\phi ,{\rm{ }}\theta ,{\rm{ }}\psi } \right]^T}\end{align}

(18) \begin{align}{\bf{u}} = {\left[ {A,\omega } \right]^T} + {\bf{w}}\;\end{align}

(19) \begin{align}{\bf{z}} = {\left[ {{U_{t,}}\alpha ,{\rm{ }}\beta ,{\rm{ }}\phi ,{\rm{ }}\theta ,{\rm{ }}\psi } \right]^T} + {\bf{v}}\end{align}

where A = [a _x , a _y , a _z ] ^T is the output of accelerometers, and $\omega$ = [p, q, r] ^T is the rate gyroscope output of the IRU or IMU. The air data variables, including true airspeed, angle-of-attack ( $\alpha$ ), angle-of-side slip ( $\beta$ ) and attitude angles ( $\phi$ , $\theta$ , $\psi$ ;) are represented by the state vector (x). The input vector (u) represents the IRU output along the body axes. The measurement vector (z) represents the measured velocities along the x, y, and z axes transformed into body axes and attitude angles provided by the INS/GPS navigation system. The process noise vector (w) consists of three elements, such as the input noise components of accelerometers (w ^A ), input noise components from rate gyros (w $^{\omega}$ ) and the random process noise of wind velocities (w ^Vw ), which is given by

(20) \begin{align}{\bf{w}} = {\left[ {{{\bf{w}}^A} \quad {{\bf{w}}^\omega }\quad {{\bf{w}}^{Vw}}} \right]^T}\end{align}

Due to the random nature of the unknown disturbance (turbulence), it is difficult to predict the local wind field variations. Moreover, considering wind velocity as a state will increase the computational burden. The wind velocity due to turbulence is considered as a random walk process model. Note that, the wind velocity due to turbulence is considered as an unknown input to the system, and it can be isolated when there is no requirement of existence conditions. It is assumed to be a Gaussian random process whose covariance can be estimated. When the noise level of the process varies intensely, an optimal wind estimation is required to calculate accurate air data. Hence, the random variation of wind velocities due to turbulence was captured using EW-AEKF and is represented as

(21) \begin{align}{\dot V_w} = {[{\dot V_{wx}}\,{\dot V_{wy}}\,{\dot V_{wz}}]^T} = {{\bf{w}}^{Vw}}\end{align}

where V _wx , V _wy , and V _wz are the components of wind velocities represented in the earth frame, and they can be transformed into the body frame using the transformation matrix. The proposed work assumed that all the sensors were fault-free throughout the PADS-FDD process, except the redundant air data sensors. Specifically, the operation of the inertial measurement unit (IMU) that includes the rate gyro sensors (p, q, r) and accelerometers (a _x , a _y , a _z ) is assumed fault-free throughout the PADS-FDD, and the noise variations of all the sensors are reduced by the EW-AEKF. In wind disturbance conditions, the EW-AEKF adaptively tuned the noise covariance of the IMU. Hence, the noise effects of rate gyros and accelerometers are well-tackled for the reliable estimation of air data, including U _t , $\alpha$ and $\beta$ .

3.1 Estimation framework of EW-AEKF

The optimality of traditional extended Kalman filter (EKF) algorithms depends on the a priori noise statistics of process noises and measurement noises (Q and R) of sensors. Prabhu and Anitha [Reference Prabhu and Anitha26] presented the air data FD using an iterated optimal extended Kalman filter (IOEKF) with the standard deviations provided by the manufacturer of the sensors. However, in the case of large variations in the sensor measurements caused by external disturbances, the known process noise (Q) cannot provide satisfactory performance for the air data estimation using IOEKF. Usually, the choice of Q _k-1 made by the trial-and-error method may induce additional uncertainty due to the user’s expertise and background [Reference Mukherjee, Das and Sadhu27]. As such, the selection of the process noise covariance matrix is a challenging part of using the traditional EKF algorithms. To resolve this issue, an adaptive version of EKF is proposed that can automatically tune the noise statistics of sensors. In this case, a computationally effective nonlinear filter such as the exponentially weighted adaptive extended Kalman filter (EW-AEKF) is considered for the estimation of air data to present the PADS-FDD algorithm. The more general form of an EKF filtering algorithm is presented as follows:

1. Step 1: Initialisation

   (22) \begin{align}\left\{ \begin{array}{l}{{{\hat{\textrm{x}}}}_0} = E({{\rm{x}}_0})\\ \\[-8pt]{{\rm{P}}_0} = E(({{\rm{x}}_0} - {{{\hat{\textrm{x}}}}_0}){\rm{ (}}{{\rm{x}}_0} - {{{\hat{\textrm{x}}}}_0}{)^T})\end{array} \right.\end{align}

1. Step 2: Prediction

   (23) \begin{align}\begin{array}{l}{{{\hat{\textrm{x}}}}_{k/k - 1}} = {{\rm{F}}_{k,k - 1}}{{{\hat{\textrm{x}}}}_{k - 1}} + {{\rm{G}}_{k - 1}}{{\rm{u}}_{k - 1}}\\ \\[-8pt]{\rm{P}} = {{\rm{F}}_{k,k - 1}}{{\rm{P}}_{k - 1}}{\rm{F}}_{k,k - 1}^T + {{\rm{Q}}_{k - 1}}\end{array}\end{align}

1. Step 3: Correction

   (24) \begin{align}\begin{array}{l}{{\rm{K}}_k} = {{\rm{P}}_{k/k - 1}}{\rm{H}}_k^T{[{{\rm{S}}_k}]^{ - 1}}\\ \\[-8pt]{{{\hat{\textrm{x}}}}_k} = {{{\hat{\textrm{x}}}}_{k/k - 1}} + {{\rm{K}}_k}[{{\rm{v}}_k}]\\ \\[-8pt]{{\rm{P}}_k} = ({\rm{I}} - {{\rm{K}}_k}{{\rm{H}}_k}){{\rm{P}}_{k/k - 1}}\end{array}\end{align}

In Equation (22), ${\hat{\bf{X}}_0}$ and ${{\bf{P}}_0}$ are the estimated initial state and state error covariance matrices, respectively. In Equation (23), $\;{{\bf{F}}_k}$ is the state transition matrix, ${{\bf{G}}_k}$ is the control input matrix; ${{\bf{u}}_k}\;$ is the control input vector; and ${{\bf{Q}}_k}$ represents the process noise covariance matrix. In Equation (24), ${{\bf{K}}_k}$ is the Kalman gain matrix, ${{\bf{H}}_k}$ is the measuring matrix. ${{\bf{S}}_k}$ is the innovation covariance, and ${{\bf{v}}_k}$ represents the measurement noise with the covariance matrix, R _k . The residual ${e_k}$ is defined as the difference between true measurements and its estimated values available at time step k, and it can be presented as

(25) \begin{align}{e_k} = {{\rm{z}}_k} - {{\rm{H}}_k}{{\hat{\textrm{x}}}}_{{{k/k - 1}}}\end{align}

The theoretical covariance of innovation sequence (C _ek ) is presented by

(26) \begin{align}{{\rm{C}}_{ek}} = E\left[ {{e_k}e_{_k}^T} \right]\end{align}

Using the maximum likelihood estimation (MLE) approach, the noise covariance can be estimated and adjusted adaptively from the measurements with the maximum probability approach. As stated in the MLE procedure, if k < N, the ideal value of the estimated ${{\bf{C}}_{ek}}$ in a moving time window is presented by

(27) \begin{align}{\hat{\textrm{C}}_{ek}} = \frac{1}{N}\sum\limits_{j = {j_0}}^k {{e_j}e_j^T} \end{align}

where N is the size of the moving time window, ${j_0}$ = k-N+1 is the first stage inside estimation, and j is the moving counter inside the estimation window. In general, the determination of system noise aims to guarantee the accuracy of the filtering process. To overcome estimation error due to external disturbances (turbulence), a moving time window (N) for robust and adaptive estimation of air data that is insensitive to variation in the monitoring is process used. Note that the proposed method requires tuning a Q matrix whose noise statistics are unknown. From the proposed system knowledge, it is clear that there are possibilities for unexpected external disturbances to the system, such as atmospheric turbulence. Therefore, a moving window with defined length N is considered. If k > N, the newest information will be weighted to locate the significance of newly received data. Furthermore, the weight value will be attenuated exponentially. Hence, the improved version of the ${\hat{\bf{C}}_{ek}}$ is written as,

(28) \begin{align}{\hat{\textrm{C}}_{ek}} = \lambda {\hat{\textrm{C}}_{ek - 1}} + \frac{{1 - \lambda }}{{1 - {\lambda ^N}}}({e_k}e_{_k}^T - {\lambda ^N}{e_{k - N}}e_{_{k - N}}^T)\end{align}

where λ is the attenuation factor. The improved ${\hat{\bf{C}}_{ek}}$ will be placed into the Kalman gain matrix to strengthen the robustness of the EW-AEKF because of the filter gain tuned to the noise fluctuations. Hence, the adaptive estimate of the process noise covariance matrix is given by

(29) \begin{align}{\hat{\textrm{Q}}_k} = {{\rm{K}}_k}{{\rm{C}}_{ek}}{\rm{K}}_{_k}^T\end{align}

Based on the adaptive Q matrix estimation, the filter gain (K _k ) was improved to strengthen the robustness of the filtering algorithm. Note that the proposed method considers tuning of the Q matrix only. Tuning of the R matrix is not required for the proposed algorithm as the measurement sensor’s noise covariance is known, reliable, and less affected by external disturbances. Furthermore, avoiding the R matrix tuning excludes the additional computations.

4.0 PADS-FDD algorithm

In this section, a simple fault detection (FDD) logic is presented to detect any type of fault on the primary air data sensors (PADS) considered for the proposed work. The FDD scheme is shown in Fig. 2.

Figure 2. PADS-FDD scheme.

In general, residual-based FDD carries the fault information irrespective of fault types. In this study, the residuals for the FDD scheme are calculated offline to simulate various types of fault scenarios by injecting faults into the fault-free flight data measurements. In order to detect a fault, the residuals are used as a fault detector or vector. The residuals calculated for the proposed PADS-FDD algorithm were calculated using the actual or redundant physical sensor measurements and the estimated data, i.e. the analytically redundant measurements. For better FDD performance, a fault model is framed, which is a formal representation of the knowledge of possible faults and how they influence the process. As far as air data sensor fault detection is concerned, the fault types will vary depending on the various flight conditions and sensor positioning on the aircraft. For example, Pitot tube ports are usually blocked by foreign objects such as the entry of insects when aircraft are on the ground, icing at high altitudes, volcanic ashes, rainwater entry, etc. Note that in all these cases, the Pitot tubes will be blocked, which may cause frozen output or a stuck fault to occur. A fault model (f _k ) includes each sensor’s data for fault types such as stuck, bias, drift and increased noise faults under abnormal conditions, including adverse weather conditions and atmospheric turbulence. The faulty sensor output can be presented as

(30) \begin{align}{Z_k} = {z_k} + {f_k}\end{align}

where Z _k is the faulty sensor data, with an added fault model at time t being the fault-free sensor data ( ${z_k}$ ). It is assumed that the estimated air data is free of faults throughout the FDD process. A residual is a computed quantity that is used to alarm if a fault is present in the supervised air data sensors (PADS) or not. The residuals (e _k ) of each sensor category of the air data system are calculated using the more general form of the residual equation as follows:

(31) \begin{align}\left| {{e_k}} \right| = \left| {{Z_k} - {{\hat z}_k}} \right|\end{align}

where $\left| {{e_k}} \right|$ is the absolute value of the residual signal, which can be explored using Equation (2) to differentiate each category of faulty sensors in the air data system. It can be presented as,

(32) \begin{align}\left\{ \begin{array}{l}\left| {{e_{{U_t}}}} \right| = \left| {({U_t} + {f_i}) - {{\hat U}_t}} \right|\\ \\[-8pt]\left| {{e_\alpha }} \right| = \left| {(\alpha + {f_i}) - \hat \alpha } \right|\\ \\[-8pt]\left| {{e_\beta }} \right| = \left| {(\beta + {f_i}) - \hat \beta } \right|\end{array} \right.\end{align}

where ${U_t}$ , $\alpha $ , and $\beta $ are the physical air data sensor measurements such as airspeed, AOA and AOSS, respectively, and f _i represents the fault added to the i ^th physical sensor measurements at time t. In order to detect a fault, a predetermined threshold (J _T ) is assigned to each sensor. The objective of the proposed method is to minimise the number of false alarms and avoid false detections in fault-free conditions. When setting a fault detection threshold, it is crucial to keep the false alarm rate (FAR) and the missed detection rate (MDR) under control. For example, a too-low threshold value will lead to an increased false alarm rate, and a too-high threshold will lead to more missed detections [Reference Chen and Patton28]. Note that increased false alarms indicate a poorly designed mathematical model, which may enable unwanted control actions that could adversely affect the system’s performance. The rationale behind this philosophy is that we only focus on detecting faults that may help to design a reliable FDD for the application of aircraft control systems related to the air data system. An alarm signal (F_A) is used to indicate faulty situations using binary decision logic (0 and 1). When a fault occurs, the F_A will set logic 1. When there are no faults, the F_A will be set to logic 0. If any of the air data sensor residual values exceed the J _T , F_A will be set to 1, and a fault is detected. It can be represented as

(33) \begin{align}&\left| {e_k^i} \right| \lt {J_T},\,{\rm{Fault}}\,{\rm{alarm, }}\ {{\rm{F}}_{\rm{A}}}\ {\rm{set}}\ {\rm{to}}\,0 \nonumber\\[3pt]&\left| {{\rm{e}}_{\rm{k}}^{\rm{i}}} \right|{\rm{ \gt }}{{\rm{J}}_{\rm{T}}}{\rm{,}}\,{\rm{Fault}}\,{\rm{alarm,}}\,{{\rm{F}}_{\rm{A}}}\,{\rm{set}}\,{\rm{to}}\ 1\end{align}

where $\left| {e_k^i} \right|$ is the absolute (modulus) of the residuals of i ^th sensor of the air data cluster (U _t , $\alpha$ , $\beta $ ) considered to be a non-negative quantity without regard to its sign. In Section 5, we will present the various fault scenarios of the sensors to simulate the proposed PADS-FDD algorithm.

5.0 Experimental setup and simulation results

In this study, we used the X-Plane 10 flight simulator to carry out an experimental flight and generate flight sensor data for the proposed PADS-FDD. The simulation setup consists of the X-Plane version 10.32 simulator software, a Logitech Extreme Pro Joystick, a Mobile Workstation Laptop with a Core i7-6820HQ CPU, and an Nvidia Quadro graphics card used as shown in Fig. 3. X-Plane is a comprehensive and powerful flight simulator that can simulate the flying characteristics of aircraft accurately. In addition to being used for certified flight simulations in accordance with FAA regulations, X-Plane has been the subject of numerous studies [Reference Bittar, Figuereido, Guimaraes and Mendes29–Reference Silva, Da Silva and Grundling32]. The turbulence generated by the X-Plane simulator is based on overlaid sinusoidal waves with various frequencies. The aircraft models provided in the X-Plane simulator are based on the blade element theory (the airplane is divided into many small elements, and the forces acting on each element are calculated with a frequency of several hertz). It includes the effects of atmospheric disturbances such as turbulence, wind shear, wind gust, etc. For different layers of the atmosphere, we can set the wind speed, its direction and the turbulence level. X-Plane can transmit and receive flight data using the User Datagram Protocol (UDP). The data I/O window of the X-Plane is employed to choose the parameters that will be communicated to Simulink in the form of data packets. UDP is used for data transmission between Simulink and X-Plane. We ran the simulator via UDP communication incorporated in Matlab/Simulink and recorded flight data under various turbulence conditions. To obtain flight data, a software-in-the-loop (SITL) simulation was performed with a Boeing B737-800 aircraft that is readily available in the X-Plane’s aircraft folder. The B737-800 aircraft is a large fixed-wing commercial transport aircraft with twin jet engines. To complete the autopilot’s fundamental operational loops, including pitch control loop (PCL), roll control loop (RCL), altitude hold (AH) and heading hold (HH), the controller receives feedback from the X-Plane flight data. The PID controller was used to implement each of these loops. UDP send and receive blocks were accessed with two independent unidirectional transmissions and an Internet Protocol (IP) address. For the UDP connection between Matlab/Simulink and X-Plane, the IP address used was 127.0.0.1 with the output ports 49,000 and 49,004. The frame rate was set to 50 packets per second (pps) during flight simulation.

Figure 3. X-Plane flight simulator setup.

The flight data was collected with the required variables. Data outputs are recorded by using the data output option in X-Plane with a text (.txt) file. The flight simulation was carried out with take-off, climb and cruise phases under atmospheric turbulence. The atmospheric turbulence magnitude was set to 2 using the weather setting under the environment tab, as shown in Fig. 4. Before implementing the proposed approach, the accuracy and stability of the filtering mechanisms were studied to verify the air data estimation algorithms using traditional EKF and the proposed EW-AEKF.

Figure 4. Environmental setup using X-Plane simulator.

5.1 Filter estimation performance results

In this subsection, we present the simulation test results of conventional EKF and adaptive EKF (EW-AEKF) to compare and verify the estimation performance. For the estimation of air data, the filtering performance of a so-called non-linear filter (EKF) and EW-AEKF was tested under abrupt environmental variations such as atmospheric turbulence. The EKF performance was found to be insufficient due to convergence errors. Therefore, an adaptive version of EKF (EW-AEKF) is preferred to estimate air data and detect sensor faults. Figures 5 and 6 show the estimation plots of air data (true airspeed, angle-of-attack, and angle-of-sideslip) under fault-free conditions. To test the filters, the simulated flight data was logged with moderate atmospheric turbulence (turbulence magnitude of 2) at an altitude of 5,000m.

Figure 5. Air data estimation using EKF.

Figure 6. Air data estimation using EW-AEKF.

As can be seen, the estimation performance of EW-AEKF is better than the estimation using EKF in terms of filter convergence error. Note that the filter estimation error is crucial for model-based FDD methods to ensure reliable estimates of the aircraft parameters. The estimation performance was compared with the root mean square error (RMSE) of the filters and is presented in Table 1.

Table 1. Comparison of filter estimation performance

5.1.1 Robustness analysis using Monte Carlo simulations

In this section, experiment results using Monte Carlo Simulations (MCS) are presented for the performance evaluation of the developed FDD scheme with respect to the uncertainty acting on the system. The analysis using MCS represents a further method for estimating the reliability and robustness of the FDD scheme [Reference Benini, Bonfù, Castaldi, Geri and Simani33]. The simulation of different fault-free and faulty data sequences was carried out with various turbulence conditions. The MCS was performed by varying the parameters of the kinematic model developed for the proposed algorithm by exploiting the MATLAB platform. Table 2 summarises the results obtained by considering the EW-AEKF estimation for the air data sensor FDD. As shown in Fig. 7, the mean and standard deviations of estimation error are within the acceptable level. The performance of the FDD method was evaluated on a number of Monte Carlo runs equal to 2,000. Based on the outcome of MCS, performance metrics such as false-alarm rate (FAR) and missed detections are analysed for the fault-free and each fault case. With the proper design of the fault detection filter and choice of the FDD thresholds, we achieved FDD with better detectability and a FAR of 0.8% for airspeed, 0.5% for AOA, and 1.75% for AOSS. The results demonstrate that the robustness of the proposed FDD method can be effectively analysed using Monte Carlo simulations.

Table 2. Monte-Carlo outcome of air data estimation and FDD performance

Figure 7. MCS results of estimation error and FAR.

5.2 PADS-FDD results using simulated flight sensor data

In this subsection, we presented the results of the PADS-FDD algorithm using simulated flight sensor data of a B737-800 aircraft with a turbulence magnitude of 2 at an altitude of 10,668m. To validate the algorithm using captured flight data, the EW-AEKF estimated air data, and the physical sensor data of the aircraft used. Note that the actual air data sensor measurements are not used in the filtering procedure and are used only as a comparison reference with the estimated air data to detect and diagnose faults. The PADS-FDD results for the various fault scenarios under atmospheric turbulence were considered. The most common fault that occurs on Pitot tubes is the blockage fault, which leads to a stuck or frozen output on the airspeed indications. Figures 8 and 9 show the Pitot tube blockage (stuck) fault in the flight sensor data at time t = 66sec.

Figure 8. Airspeed sensor stuck fault at t = 66sec.

Figure 9. Airspeed sensor residual during stuck fault.

As can be noticed, the absolute value of the residual signal ( ${e_U}_{_t}$ ) exceeds the predetermined threshold of 2m/s at time t = 66sec, and it jumps beyond 30m/s. Hence, the fault is detected. During the sensor fault, the RMSE of the filter (EW-AEKF) is 15.7612m/s. Note that the fixed threshold value is related to the permissible error limit of an airspeed indicating system as mentioned in the EASA specifications. Figures 10 and 11 show the simulation results of the drift fault and residuals of the airspeed sensor.

Figure 10. Airspeed sensor drift fault at t = 70sec.

Figure 11. Airspeed sensor residual during drift fault.

In Fig. 10, the drift fault of the airspeed sensor is presented by injecting the fault in the physical airspeed sensor measurements at time t = 70sec.

Drift fault is a condition where the measurement error rises over time, which originates from the loss of the sensor’s sensitivity. In Fig. 11, the residual ( ${e_U}_{_t}$ ) exceeds the 2m/s threshold and reaches beyond 3.5 m/s, which is not acceptable. The RMSE of the filter during the airspeed sensor drift fault is 2.5129m/s. In Figs 12, 13, 14, and 15, the airflow angles ( $\alpha$ and $\beta$ ) fault detection is presented as a bias fault, which is the most common fault experienced by the commercial aircraft vane type $\alpha$ and $\beta$ sensors. A bias fault is the steady offset error between the actual and measured values. The results are shown below.

Figure 12. Angle-of-attack ( $\alpha$ ) sensor bias fault at t = 65sec.

Figure 13. $\alpha$ sensor residual during bias fault.

Figure 14. $\beta$ sensor bias fault at t = 60sec.

Figure 15. Angle-of-sideslip ( $\beta$ ) sensor residual at t = 60sec.

In Fig. 12, the bias fault is introduced to the physical angle-of-attack sensor measurements at time t = 65sec. The absolute value of the residual signal exceeds the predetermined threshold of 0.008rad. As a result, the fault is detected at t = 65sec. Figure 13 shows the residual variations during bias fault on the angle-of-attack ( $\alpha$ ) sensor. Figure 14 shows the bias fault on the angle-of-sideslip ( $\beta$ ) sensor at time t = 60sec.

In Fig. 15, the residuals after the bias fault on the $\beta$ sensor, which exceed the threshold value of 0.0025rad, is presented. The RMSE of the filter during $\alpha$ and $\beta$ sensors bias faults is 0.0074rad and 0.0012rad, respectively.

In Figs 16, 17, 18, 19, 20, and 21, an increased noise fault on the air data sensors is introduced under the atmospheric turbulence. During turbulence, the air data sensors’ noise increases significantly because of the heavy airflow distortion. The turbulence magnitude of 3 was set in the flight simulator weather setup window. The results are shown below:

Figure 16. Airspeed sensor increased noise fault at t = 141.9sec.

Figure 17. Airspeed sensor residuals after noise fault.

Figure 18. $\alpha$ sensor increased noise fault at t = 141.95sec.

Figure 19. $\alpha$ sensor residuals after noise fault.

Figure 20. $\beta$ sensor increased noise fault at t = 141.9sec.

Figure 21. $\beta$ sensor residuals after noise fault.

Figure 16 shows the increased noise fault on the airspeed sensor under external disturbance (turbulence). The residual variations found in Fig. 17 show that the fault detected at time t = 141.9sec, and the residuals surpass the fixed threshold value of 2m/s. As can be noticed, the airspeed sensor residual jumps beyond 16m/s after the occurrence of the fault. The RMSE of the filter during the fault is 3.5572m/s.

As seen in Figs 18 and 20, increased noise faults occurred on $\alpha$ and $\beta$ sensor measurements at time t = 141.95 and t = 141.9sec, respectively. As can be noticed, heavy airflow distortions caused by the increased turbulence magnitude on the rotating vane type $\alpha$ and $\beta$ sensors cause an increased noise fault. In Figs 19 and 21, the AOA ( $\alpha$ ) and AOSS ( $\beta$ ) sensor residuals indicate the occurrence of increased noise faults when they exceed the fixed thresholds of 0.01 and 0.005rad, respectively. The RMSEs of the filter during $\alpha$ and $\beta$ sensor noise faults are 0.0039 and 0.0021rad, respectively.