Abstract
A Fault Detection and Isolation (FDI) algorithm is developed to protect against Water-Blockage (WB) pitot tube failure in the safety-critical Air Data System (ADS) used on small Unmanned Aircraft Systems (UAS). The algorithm utilizes two identical Synthetic Air Data Systems (SADS) as the basis for state estimation. Each SADS works independently with a pitot tube while sharing an IMU and GNSS receiver. The fault detection is designed using the integrity monitoring framework, and the isolation is obtained via independent fault detection channels. The ADS requirements are established, and the WB failure mode is analyzed based on real faulty air data. A new residual-based test statistic is introduced, and the link among the test statistic, observability matrix, and Minimal Detectable Error (MDE) are examined. Finally, a flight data set with a known water-blockage fault signature is used to assess the algorithm’s performance in terms of the air data protection levels and alert limits.
1 INTRODUCTION
A reliable Air Data System (ADS) plays a vital role in an aircraft’s safety and performance. While ADS provides the measurement of various parameters, airspeed Va, angle-of-attack α, and angle-of-sideslip β are the main parameters that define the flight envelope. As shown in Figure 1, airspeed is the speed of an aircraft relative to the air, and angle-of-attack and angle-of-sideslip are the flow angles relative to the aircraft. Air data is usually measured onboard by accurate air data sensors such as the pitot-static tube and the angle vane. Also, reliability analyses such as Fault Tree Analysis (FTA) and Failure Mode and Effect Analysis (FMEA) are often used to help identify fault modes and certify the redundant hardware systems.
Illustration of air data triplet: airspeed Va, angle-of-attack α, and angle-of-sideslip β
The analytical redundancy approach is particularly useful for small Unmanned Aircraft Systems (UAS) due to the Size, Weight, And Power (SWAP) requirements. However, there is a lack of analytical methods to certify analytical redundancy. As emerging technologies such as Urban Air Mobility (UAM) (Vascik et al., 2018), or UAS operations either in the Line-of-Sight (LOS) or Beyond Visual Line-Of-Sight (BVLOS) (Cour-Harbo, 2017; Fang et al., 2018; Johnson et al., 2017; McCrink & Gregory, 2018; Yapp et al., 2018) mature, the need for rigorous and certifiable analytical redundancy methods will increase. In these applications, ADS is one of the safety-critical subsystems which needs to be certifiable and meet safety requirements.
Most of small UAS usually have one or two pitot tubes as the only sensors in its ADS. For example, the recent Part-135 certification process requires any unmanned air carriers (i.e., package delivery) to have at least one heated pitot tube (Federal Aviation Administration, 2019). However, many small UAS cannot afford to have a heated pitot tube onboard due to its cost and power requirements. Without the heating system, the low-cost pitot tubes on many small Unmanned Aerial Vehicles (UAVs) are prone to Water-Blockage (WB) faults. This is why many small UAV operations, such as agricultural surveying and construction inspection, cannot be carried out reliably during the rainy days. In Figure 2, a typical inexpensive pitot tube [10 US dollars to 20 US dollars (JDrones, 2020; Eagle Tree Systems, 2020)] is shown. It can be seen that the pitot tube is connected to a transducer via plastic tubes. The setup is simple and used by many UAVs but prone to the WB faults since there is no built-in drainage or heating system in the pitot tube. Water can enter the pitot tube on flights during foggy or rainy days, which fully or partially block the stagnation ports and affect the transducer’s pressure readings.
Entire ADS using a pitot-static tube, pressure tubes, and a pressure transducer for a small UAS
To improve the safety and reliability of the ADS and minimize the number of redundant and multiple sensors, one approach being considered is called a Synthetic Air Data System (SADS) (Lie & Gebre-Egziabher, 2013; Sun et al., 2019b). A SADS is an estimator that calculates air data quantities using non-air data sensors such as the GNSS, IMU, magnetometer, and mathematical model of the aircraft.
SADS is a form of analytical redundancy that can help detect and deal with faults in the traditional ADS of small UAS. SADS can also potentially be coupled with one or two air data sensors to resolve the low-reliability issue. In fact, SADS has already been implemented in some commercial aircraft such as the Boeing 787 (Austrailian Transport Safety Bureau, 2015). The use of SADS is also being considered by many other aircraft designs at this time (Federal Democratic Republic of Ethiopia, Ministry of Transport, Aircraft Accident Investigation Bureau, 2020; Komite Nasional Keselamatan Transportasi, Republic of Indonesia, 2018; SeekingAlpha, 2019).
In what follows, we give a brief overview of the prior work on air data Fault Detection and Isolation (FDI) in the literature. We also explain why the existing methods are not adequate to certify ADS on small UAS and why the Integrity Monitoring (IM) framework can potentially solve this problem.
1.1 Prior work
Air data FDI using advanced control and estimation algorithms has renewed interest over the last decade due to the recent advancements in the safety-critical UAV applications. These air data FDI techniques can be roughly separated into three categories: model-based, model-free, and data-driven algorithms.
The model-based algorithms typically leverage the dynamic model of the aircraft (Freeman et al., 2013; Hansen & Blanke, 2014; Ossmann et al., 2017). For example, Freeman et al. (2013) designed an airspeed fault detection algorithm using the aerodynamic model of the aircraft as well as linear robust H∞ filters to detect faults, reject disturbance, and provide robustness to the modeling errors.
The model-free algorithms mainly rely on the sensor information and the kinematic models of the vehicles (Eubank et al., 2010; Guo et al., 2018; Lu et al., 2016; Van Eykeren & Chu, 2014). The model-free algorithms also often use Kalman Filter (KF)-based estimation techniques and the innovation χ2 test to determine faults. An illustrative example is shown in Lu et al. (2016). They use an adaptive three-step unscented KF to detect and isolate air data faults.
The data-driven algorithms primarily rely on large amounts of data to develop reliable input-output methods (e.g., autoregressive models, neural networks). These models are then used to detect faults through inconsistency check (Borup, 2018; Fravolini et al., 2017, 2018; Rohloff et al., 1999).
These air data FDI algorithms are primarily concerned with how to detect faults accurately, however, they do not provide a framework for ensuring the reliability of the air data fault detection algorithms from a requirement point of view. That is, they do not address the following question: Can we design an air data FDI algorithm to satisfy a given set of system requirements such as integrity and continuity, and provide statistical protection levels for the air data estimates?
Recent work (Freeman, 2014; Hu & Seiler, 2015; Kotikalpudi et al., 2020) has made progress towards certification of analytically redundant systems via reliability analysis. In the work here, we borrow tools often used in the field of integrated GNSS navigation to help design air data fault detection algorithms to ensure reliability.
One of the standard techniques for fault detection of the safety-critical aerospace navigation systems is called Receiver Autonomous Integrity Monitoring (RAIM) (Brown & Chin, 1998; Lee, 1986; Parkinson & Axelrad, 1988; Sturza, 1988). RAIM methods are used for safety-critical applications such as GNSS-based precision landing systems for aircraft (Khanafseh et al., 2014; Tanιl et al., 2017a, 2017b; Walter et al., 2008) and have been the subject of a significant amount of work for the last two decades. RAIM methods have also been used recently to provide Integrity Monitoring (IM) for other navigation systems such as the Simultaneous Localization and Mapping (SLAM) problem (Arana et al., 2019a, 2019b; Bhamidipati & Gao, 2019).
The basic idea of RAIM is to leverage redundant measurements at every time step [see snapshot detection scheme (Brown & Chin, 1998; Parkinson & Axelrad, 1988; Sturza, 1988)] or sequentially (Joerger & Pervan, 2013) to come up with probabilistic measures to detect faults and provide statistical bound to protect the state estimate. The advantage of this approach is that it provides the means for rigorous integrity risk computation. It uses redundant measurements to achieve fault detection capability and quantify the impact of undetected faults on state estimation errors. Another advantage is that the calculation of the threshold is based on probability, not selective tuning. Unlike RAIM, in many of the aforementioned works, it is often seen that a particular threshold is handpicked for a given application without rigorous probabilistic calculation.
However, the rigorous IM framework has not always been implemented on emerging non-PNT applications. This is partly due to the sensor measurements’ inhomogeneity or the non-linearity of dynamics in many applications such as UAS. Many systems, such as ADS, have limited redundant and heterogeneous measurements at every time step. This limitation sometimes makes the snapshot of residual-based detection function infeasible to determine faults. And while linearization errors are usually small in the GNSS applications, measurement models in other systems are generally highly nonlinear, and the linearization errors’ can be large. Therefore, the non-linearity might have a significant effect on the existing IM techniques. Lastly, many systems have observability issues [unobservable states (Kassas & Humphreys, 2014) or conditionally observable states (Sun et al., 2019b)], and this is usually overlooked when dealing with IM in GNSS applications.
Another aspect of fault detection in IM is the selection of an appropriate fault detector or test statistic. The goal of a fault detector is to detect fault quickly without raising too many false alarms. For real-time applications, online fault detectors are preferred. There are many online fault detectors such as the simple residual thresholding, KF innovation χ2 test, least-squares residual-based χ2 test [or commonly referred as the snapshot RAIM (Brown & Chin, 1998; Parkinson & Axelrad, 1988; Sturza, 1988) in navigation literature], Sequential Probability Ratio (SPRT), Cumulative Sum (CUSUM) test, and Generalized Likelihood Ratio (GLR) test.
Residual thresholding is the most straightforward test as it only requires a threshold to determine whether the data has exceeded the nominal level. The KF innovation χ2 test is suitable for any KF-based state estimation, and the least-squares residual-based χ2 test is a good choice when redundant measurements are available. The three methods mentioned above deal with linear or quadratic functions of the residual.
On the other hand, both SPRT and CUSUM tests are well-known for their nonlinear stopping rules (Gustafsson, 2000). For example, the standard one-sided SPRT test requires three tuning parameters: drift v, threshold h, and reset level a. The basic idea of a one-sided SPRT test is to test whether the test statistics have drifted away significantly from the threshold. The drift parameter v is used to subtract from the test statistic to control the drift’s level, and the parameter a is used to prevent negative drift.
Similarly, the one-sided CUSUM test is the same as the one-sided SPRT test with the reset level a = 0. Several variations of both SPRT and CUSUM tests can be found in the literature. However, both tests require hand-tuning for the desired outcome. The GLR is also a powerful nonlinear test for fault detection, but it usually requires the knowledge of Probability Density Functions (PDF) under different hypotheses. In this paper, the KF-based detectors will be utilized with some additional novel improvement.
1.2 Contribution
This paper provides four main contributions to the air data FDI literature: First, we design a dual pitot tube air data fault detection and isolation system that can be easily implemented on most UAVs. Second, we expand on sequential IM techniques in the Kalman filter setting to evaluate the integrity risk for the designed fault detection algorithm.
Specifically, we show how to deal with the limited redundant measurement problem and establish an analytical relationship among the residual-based test statistic, the Linear Time-Varying (LTV) observability matrix, and the Minimum Detectable Error (MDE). We also show how monitoring the observability of the system can potentially help rule out false alarms. Furthermore, we generalize the IM performance trade-off design procedure so that we can use it to evaluate other pitot tube failure modes.
Third, we also show how to establish alert limits and protection level bounds for the angle-of-attack and sideslip states. Lastly, we demonstrate our algorithm’s capability using a recorded flight data in which a known WB pitot tube failure occurred.
1.3 Paper organization
The remainder of this paper is organized as follows. Section 2 presents a brief description of the model-free SADS estimator used for the dual pitot tube air data system design developed in this paper. Section 3 presents the air data system requirements needed for the fault algorithm design. Section 4 describes the WB failure mode used in this work. Section 5 presents the fault detection design and analysis, which includes the derivation of the residual-based test statistic and its relation to the observability matrix, the MDE design and analysis, and the IM performance trade-off design procedure. Section 6 derives the alert limits and protection level calculations of angle-of-attack α and sideslip β. Section 7 presents the flight results and its associated detection performance. Concluding remarks and future outlook are given in Section 8.
2 DUAL PITOT TUBE AIR DATA SYSTEM DESIGN
2.1 System architecture
For the development that follows, we propose a dual pitot tube air data system for small UAS. The architecture consists of two identical SADS estimates of α and β by fusing airspeed measurements from the pitot tube with information from an IMU and a GNSS receiver (Sun et al., 2018, 2019a, 2019b).
Each SADS utilizes its own pitot tube, but both SADS share a GNSS receiver and IMU. Sharing the GNSS receiver and IMU reduces cost and software complexity. The design can be easily expanded to architecture with dual GNSS receivers and IMU units.
Each SADS can detect faults independently (i.e., identify and isolate the fault source) and is designed to satisfy the given system performance requirements (i.e., integrity and continuity requirements). The two SADS together provide recovery capability for a single faulty pitot tube failure via simple decision logic. This ADS system also provides accurate estimates and protection levels for the synthetic angle-of-attack α and sideslip β. The entire dual ADS design is illustrated graphically in Figure 3.
Dual air data fault detection design
In comparison to the commonly used triple-redundancy ADS design (Yeh, 1996), the dual pitot tube fault-tolerant ADS has a unique advantage: it only includes two small and inexpensive pitot tubes by leveraging a so-called dynamic redundancy approach (Isermann, 2005), which can be easily installed on the UAVs. In the case of a single pitot tube fault, the dual ADS can shut off the faulty pitot tube and continue its nominal operation using the secondary pitot tube (sometimes referred to as a hot standby).
Additionally, the two independent SADS filters can be implemented asynchronously on the hardware. It is a more fault-tolerant design choice in comparison to other filter methods [e.g., solution separation method (Joerger et al., 2014)], which in this case would use both pitot tube measurements simultaneously. Note that failure of a single pitot tube would lead to potential loss of control even if multiple filters (e.g., a bank of KF filters for various pitot tube failure modes) were used. Though simultaneous failure of both pitot tubes could occur under the same rainy condition, we considered that simultaneous-failure case beyond the scope of this paper and can be considered in future work. Table 1 summarizes the decision logic for the dual pitot tube ADS fault detection design.
Dual Pitot tube ADS decision logic for each scenario
Since SADS is the one of key components in the proposed algorithm, we will briefly go over some details to facilitate the understanding in the next subsection. For more details, please refer to Sun et al. (2019b).
2.2 Synthetic Air Data System (SADS)
The synthetic SADS estimator is an extension of the 15-state, loosely-coupled INS/GNSS EKF (Gleason & Gebre-Egizabher, 2009), which blends information from an IMU and GNSS receiver. The INS/GNSS filter’s state vector is augmented by three additional states representing the components of the wind velocity vector. Therefore, the SADS filter states, expressed in the error state vector δx ∈ ℝ18×1, is given by:
1
where δp = [δL δλ δh]T is the position error vector in latitude, longitude, and altitude, δvn = [δVN δVE δVD]T is velocity error vector resolved in the North-East-Down (NED) frame, denoted by the superscript n. The vector 
represents the attitude errors which are defined to be the small rotation angles between the actual NED frame and the estimated NED frame. The subscript nb indicates the positive direction is defined as being from the NED frame (n-frame) to the body frame (b-frame). The vectors δba = [δbax δbay δbaz]T and δbg = [δbgx δbgy δbgz]T are the accelerometer and rate gyro triad output bias error vectors, respectively. Finally, δWn = [δWN δWE δWD]T is the error in the wind velocity vector resolved in the NED frame.
The synthetic SADS estimator synthesizes an estimate of α and β without using the α and β sensor measurements. The synthetic estimates of α and β are calculated using the EKF state estimates as follows:
2
where:
3
The 
is the coordinate transformation from NED to the body frame. The measurement vector z ∈ ℝ7×1, shown in Equation (4), consists of position p and velocity vn estimates from the GNSS receiver, along with the true airspeed Va estimate determined using the pressure measurements from the pitot tube.
4
The time and covariance update equations for this filter are, for the most part, identical to those of the filter described in Gleason and Gebre-Egizabher (2009). What is new is the dynamic model for the augmented states (wind) and the measurement model.
Similar to the modeling of the accelerometer and gyroscope biases in the filter, the dynamics of the wind are modeled as a first-order Gauss-Markov model, motivated by Berman and Powell (1998). The details of the Gauss-Markov model for the wind and sensors can be found in (Berman & Powell, 1998) and (Xing, 2010), respectively. However, for the sake of completeness we re-state the process noise matrix Rw that accounts for accelerometer, gyroscope, and wind, respectively, for a quick reference:
5
where 
and 
are the accelerometer and gyroscope white noise variances, respectively. The parameter 
, 
, and 
are the accelerometer, gyroscope, and wind random walk variances, respectively. The τad, τgd, and τWd are the associated time constants defined in the first-order Markov process.
The linearized measurement model used by the EKF is δzk = Hkδxk + υk, where υk, the measurement noise vector, is assumed to follow a normal distribution with zero mean and covariance R, denoted as N(O, R). The measurement Jacobian Hk ∈ ℝ7×18 is given by:
6
where the first two block rows map the EKF states into the GNSS position and velocity measurement errors, and the last block row maps the EKF states into the airspeed measurement error. The matrix 
is derived from linearizing the nonlinear airspeed measurement model 
. The matrix 
is similarly derived and happens to be equal to 
. The matrix 
is shown in the following:
7
The associated measurement noise covariance R is shown as follows:
8
where the diagonal of R contains the position, velocity, and airspeed noise variances.
2.3 Observability consideration
An advantage of this SADS estimator is that it does not use the aircraft dynamic model, and it provides synthetic α and β estimates as well as their covariances. Specifically, unlike the model-based SADS, which uses the aerodynamic model of the aircraft and six degree-of-freedom dynamic equations, this mode-free SADS estimator only relies on the kinematic equation and sensor measurements. However, this estimator is conditionally observable as analyzed in detail in Sun et al.’s earlier work (2019b). Briefly, ensuring observability of this estimator requires the following two conditions (i.e., conditionally observable):
The airplane must be accelerating so that the INS/GNSS heading and gyro bias states become observable (Gleason & Gebre-Egizabher, 2009)
The wind vector Wn must be quasi-static. The term quasi-static means that the variations in Wn are assumed to be negligibly small over a small time window whose size is defined in (Sun et al., 2019b).
The second condition is required to ensure that changing airspeed and the wind states Wn separately observable (i.e., the wind triangle relationship). Furthermore, the quality of estimates, in part, depends on the degree of observability.
The degree of observability is determined quantitatively by analyzing the condition number of observability Gramian in Sun et al. (2019b). Since the synthetic estimate is conditionally observable, the ability to detect air data system faults is also conditional. One of this paper’s key contributions is to show how observability is related to the fault test statistic, which is explained and demonstrated in Sections 5 and 7, respectively.
3 AIR DATA SYSTEM REQUIREMENT
To quantify the air data fault detection performance, we start with given system requirements, such as integrity risk Ireq and continuity risk Creq. The integrity risk is the probability that a hazardous fault goes undetected, and continuity risk is the probability that an alarm is issued about the presence of a fault when in fact there is no fault. Mathematically, they are defined as (Pervan, 1996):
9
10
where P is shorthand for probability of, and MD, FA, DF, MI, and NI stand for missed detection, false alarm, detected failure, missed-identified failure, and nonisolable failure, respectively.
Since the primary focus here is to deal with fault detection against pitot tube faults only, the second terms on the right-hand side of Equations (9) and (10) are ignored. These terms are associated with the isolation problem, which deals with all possible fault modes associated with the pitot tube.
In our case, we assume the pitot tube only experiences one particular failure mode, so the fault detection and isolation are essentially achieved at the same time. We also rely on the simple decision logic (Table 1) to recover from the faulty pitot tube experiencing the WB failure. Therefore, we limit our scope to the following performance requirement:
11
12
As stated in Equations (11) and (12), the probability of missed detection PMD represents the probability of not detecting a pitot tube failure given the pitot tube has indeed failed. Similarly, the probability of false alarm PFA is the probability of issuing an alarm when there is no pitot tube failure. Also, if the pitot tube is working properly, this hypothesis is denoted as H0. If the pitot tube is not working correctly, this hypothesis is denoted as H1. The details of the particular WB pitot tube fault mode is analyzed in Section 4. Also, since currently there is no universally accepted the numerical values of Ireq and Creq for the WB pitot tube failure to our knowledge, we treat them as the trade-off variables in the Section 5.
4 WATER BLOCKAGE FAILURE MODE
As mentioned in Section 1, low-cost pitot tubes such as the one shown in Figure 2 are susceptible to a failure model called Water Blockage (WB) during foggy or rainy days. This is a failure mode where the water particles in the air can enter through the front of the pitot tube and accumulate, leading to a reduction in total pressure either slowly or abruptly, as illustrated in Figure 4.
Illustration of the pitot static tube experiencing WB fault scenario
The airspeed Va is typically calculated based on Bernoulli’s principle as follows1:
13
where Pt is the stagnation or total pressure, Ps is the static pressure, and ρ is the air density. A partially blocked pitot tube would affect ΔP, which often leads to an airspeed drop. The size of the airspeed drop can vary significantly based on how much water is clogging the pitot tube.
Figure 5 shows a time history of two different faulty airspeed data sets from two different UAVs. The first faulty airspeed data (the top figure in Figure 5) comes from an agricultural inspection experiment. The flight data was collected by Sentera LLC.
Two faulty airspeed data sets due to WB
The UAV took off around 570 secs, but the airspeed quickly decreased at 614 secs due to the WB faulty pitot tube. The other faulty airspeed data shown here2 is reproduced from Hansen & Blanke (2014). The airspeed experiences a sharp drop at 2335 secs due to the water-clogged pitot tube.
Even though two different UAVs operate at different flight conditions (i.e., the nominal airspeed from the second set is almost three times higher than the first one), both pitot tubes experience a similar pressure drop rate. Those two airspeed drops’ slope is estimated to be 2.5 m/s2 and 3 m/s2, respectively. Although this profile could be different from different pitot tubes, we will use them as the fault profile from which we need protection for the illustration in this paper.
Ideally, a larger faulty airspeed sample size would be required to represent the WB faulty pitot tube fault characteristics. However, it is difficult to obtain faulty airspeed data due to the WB failure in flight since: 1) the precise occurrence (i.e., timestamp) of the WB fault is usually unknown, and 2) faulty airspeed data is sensitive information and generally not shared in the public domain. In fact, to the best of our knowledge, this is the first paper that utilizes more than one set of faulty airspeed data due to the WB failure mode.
The airspeed measurement model under the faulty condition shown in Equation (14) is used for the SADS estimator. The WB fault mode is modeled as an unknown linear ramp fault, denoted as 
, and the nominal airspeed is calculated by taking the euclidean norm (2-norm) of the difference between the inertial velocity vn and wind vector Wn in the navigation frame. The airspeed measurement noise 
is modeled as the white Gaussian noise. The noise variance is typically unknown after the fault occurs. Here we assume the noise variance stays the same before and after the fault:
14
Also, we assume the fault component 
affects the nominal airspeed measurement continuously after water clogs in the tube. In other words, the water is assumed to stay in and continue clogging the pitot tube. A timeline for the fault scenario is described in Figure 6.
Fault scenario: the green portion of the timeline is fault-free section and the red is the faulted case
The fault vector 
starts entering at the time k and persists for future times. The test statistics developed in Section 5 uses a sliding window of size q to detect faults. This sliding window moves forward into the red region as time continues. Note that the sliding window would start right away when the measurement update of the KF filter reaches enough measurement for its detector. Hence the fault will most likely fall into the sliding window since the water blockage fault usually happens after the takeoff.
The methodology determining the minimum detectable faulty component 
for the fault detection design is presented in Section 5. The minimum detectable airspeed fault depends on various factors, such as the measurement sampling rate Ts, the sliding window size q of the fault test statistic, the integrity risk PMD, and continuity risk requirements PFA. An IM trade-off design procedure is also presented to show how various factors can affect the design choice based on a given set of requirements.
5 FAULT DETECTION DESIGN AND ANALYSIS
In this section, we first discuss the choice of Kalman filter based test statistics. Specifically, in addition to the conventional innovation-based KF test statistic, we introduce a sequential residual-based test statistic. We show how this test statistic is related to the observability matrix. Second, we discuss how the test statistics are generally designed to meet the integrity requirements.
The Minimum Detector Error (MDE) metric is used to link the integrity requirements and the specific air data system requirements. The MDE metric is determined through the Detector Operating Characteristic (DOC) curves. We also provide a general design procedure to determine the acceptable MDE and detection time τ for different failure modes. Lastly, we examine the quality of the proposed KF residual-based test statistic through the MDE analysis.
5.1 Kalman filter based test statistics
Since the fault detection design in this paper relies on the Kalman filter based estimation, both innovation and least-squares residual-based χ2 tests are considered for the fault detector design. Kalman filter based test statistics are widely used in the field of GNSS applications.
5.1.1 Innovation-based test statistic
The most popular method is called the normalized innovation squared χ2 test (Bar-Shalom et al., 2002). It uses the innovation vector and its covariance to form a test statistic, which follows a central χ2 distribution with Df = mq under the fault-free hypothesis H0 as shown below:
15
where Dγ, k is test statistic at the time step k, γk is the innovation vector calculated using 
, the matrix Sk is the innovation covariance calculated from 
, m the number of measurements at the time step k, and q is the sliding window size. The vector 
is the predicted estimate, and 
is the covariance of the KF prediction.
A suitable threshold for this test can be computed by using the inverse chi-square cdf 
given the desired probability of false alarm and the appropriate Df as follows:
16
The innovation-based test statistic’s effectiveness depends on the quality of the KF-predicted estimates (linear prediction), KF measurements (sensor quality), and length of q. However, one weakness of the innovation-based test statistic is that it cannot be analyzed easily if there is a fault. This is because all the information embedded in the innovation vector and its covariance, and the contribution from a fault cannot be parsed out analytically. In other words, it would be more beneficial if we could find a test statistic that is both analyzable and informative.
5.1.2 Residual-based test statistic
The other common fault detector for KF-based estimation is the least-squares residual-based χ2 test. However, the residual-based χ2 test is not applicable to the SADS considered here because there is no redundant airspeed measurement at every time step; past measurements will be required for the fault detection design instead of using the well-established snapshot RAIM method. To overcome this issue, we formulate a KF residual-based test in the following:
17
where n is the number of the states in the KF and rk−q:k is stacked residual vector from the past time step k−q to the current time step k, denoted as 
. Each residual rk is computed from the difference between the measurements and a posteriori estimate 
. The matrix Σ is the covariance matrix of the weighted least residual vector rk−q:k is shown in Equation (18):
18
where ⊗ is the Kronecker tensor product and Rw is the process noise matrix. The matrix 
is realized through the batch linear system realization shown in Appendix A. The sliding window residual-based test statistic is used here instead of a one-time step residual test because the number of measurements m is less than the number of states n for the SADS.
In other words, the popular snapshot RAIM method from GNSS does not work here since no redundant measurements are available at each time step. The χ2 test requires Df = mq − n > 0, therefore the threshold for the residual-based test statistic is calculated as follows:
19
The sequential residual-based χ2 fault detection test statistic has similar properties to the snapshot RAIM method. Also, we make a connection between this residual-based test statistic and the LTV observability matrix. This is done by connecting a window of measurement to a past state vector xk−q using weighted least squares. The state vector xk−q can be estimated by applying weighted linear least squares to the batch linear system shown in Equation (A2):
20
where 
signifies the posteriori estimate from the past measurement from time step k−q to the current measurement k. The matrix 
is calculated as:
21
where 
is the LTV observability matrix shown in Equation (A7). The residual vector rk−q:k can be subsequently expressed in the following:
22
Under the fault-free H0 hypothesis, we can also write rk−q:k as follows:
23
where rk−q:k follows a normal distribution N(O, Σ), and Wk−q:k and Vk−q:k are defined in Equation (A4) and (A5) respectively.
Using Equation (22), we can also write Dr,k as follows:
24
where the subscripts mq × mq and k−q:k are dropped to shorten the notation. The last equality of Equation (24) is obtained because both Σ−1 and 
are symmetric, and the matrix 
is idempotent. The formal proof of this matrix equality is given in Appendix B.
The mathematical revelation in Equation (24) shows that the KF residual-based test statistic is a function of the observability matrix. Furthermore, the matrix 
inside of 
is the discrete weighted observability Gramian or the Fisher information matrix (Bar-Shalom et al., 2002).
This test statistic has a distinct advantage: it gives users a tool to analyze how the system’s observability affects the test statistic Dr,k given a sliding window of the measurement from time step k-q to k. By analyzing the observability Gramian, we can tell how well the current system is observable. We can make useful statements between the effect of the test statistic and the motion of the vehicle (indirectly represented by the observability matrix). Furthermore, it can also be used to eliminate false alarms, as demonstrated in Section 7. Hence, we also call this test statistic the observability-based test statistic.
5.1.3 Limitation of snapshot RAIM test statistic
It is worth noting that the test statistic in Equation (24) is different from the well-known RAIM-like 
test statistic, which also follows a central χ2 with degrees-of-freedom mq − n under H0 hypothesis.
This test statistic does not account for the process noise from the time update step in the KF prediction step. It loses all the dynamic information between the KF measurement update. Again, though the snapshot 
test statistic is often used for GNSS integrity monitoring, it is inapplicable when dealing with the system considered here which does not have redundancy measurements at each time step.
5.2 Minimum Detectable Error (MDE) design
Before we proceed with determining the Minimum Detector Error (MDE), we will discuss some concepts and define some terms that will be used later. Any fault detection algorithm’s goal is to detect credible faults before they lead a hazardous situation (e.g., loss of the aircraft, collision with terrain).
For a given UAV, the stall angle of attack αstall and the minimum airspeed at which the airplane can fly Va,stall (stall speed) are synonymous. We will assume we are dealing with an electrically powered UAV, so its mass does not change during flight. Thus, the flight detection algorithm we design will have to detect faults before the estimated airspeed falls below Va,stall.
Since the UAV’s operating speed Va is generally not constant during a given flight, the allowed drop in the estimate of airspeed (before the airplane is outside of the safe-flight envelope) is not constant either. To simplify the design, we use the average operating speed 
as an approximation. We will call the difference between average operating airspeed 
and Va,stall the airspeed Allowable Error, or AE for short.
Given a WB-fault profile, we can determine the time required for the airspeed estimate to drop below the stall speed. We call this τmax, and the fault detection algorithm must detect a WB fault in a time shorter than τmax.
Finally, we call the smallest airspeed estimation error that can be detected consistently (quantified by the missed detection and false alarm rate probabilities) the MDE. The fault detection algorithm ensures that MDE< AE and raises the alarm when the detection time less than τmax after the onset of a pitot tube failure.
5.2.1 MDE and τmax determination
By examining the slope of the faulty airspeed data in Figure 5, it is determined that the minimum faulty airspeed drop rate is about 2.5 m/s2. For the particular UAV used in the flight experiment, the average operating airspeed is about 17.5 m/s, and the stall speed is about 10 m/s. The difference 7.5 m/s between 
and Va,stall translates to the maximum detection time τmax = 3s as follows:
25
If the fault detection algorithm fails to detect the fault before the fault exceeds AE, then the algorithm is ineffective against the allowable fault. Mathematically, the AE should satisfy the following:
26
where the lower bound MDE is a function of PFA and PMD. The AE = 7.5 m/s is reasonable, but a tight bound since the UAV should be able to recover even if the airspeed drops below the stall speed as long as there is a sufficient altitude for recovery. Hence, a larger AE can be found based on the average operating altitude. Nevertheless, we will use 7.5 m/s as the upper bound for AE in the following analysis.
5.2.2 MDE and sampling rate
We model the airspeed fault as a linear ramp fault using the minimum drop rate of 2.5 m/s2. The fault detector should still catch any rate that is higher 2.5 m/s2 since a larger fault would result in a quicker detection.
Since χ2 based tests are used for the fault detection, one of the requirements for χ2 test is that the degree-of-freedom (Df) has to be greater than zero. For example, the least-squares residual-based fault detection method (Brown & Chin, 1998) requires the number of the measurement m to be greater than the state n. Because we are doing sequential measurement update for the airspeed (i.e., m = 1 at each time step k) and we only assume the fault comes from the pitot tube, the minimum number of airspeed measurements needed for the residual-based χ2 test is 19 since the number of states is 18.
Therefore, we need to collect at least 19 airspeed measurements (i.e., sliding window size q = 19) to detect fault within the maximum allowable detection time. Sampling too fast would also degrade the KF’s performance because measurements closely spaced in time would cause the innovation vector to be highly correlated with itself, which violates the uncorrelatedness assumption.
Figure 7 shows the simulated airspeed linear ramp fault profiles at 1 s and 0.08 s sampling time over 3 s using the 2.5 m/s airspeed drop rate. The sampling time Ts = 1 s might be too slow for detection as the number of faulted measurement is too small for detection test statistic to respond before the time exceeds τmax.
Simulated airspeed linear ramp fault profile at two different sampling rates
On the other hand, the sampling time Ts = 0.08s is sufficient even if we only use the measurements after the fault occurs to form the test statistic3. Therefore, an appropriate sampling time should be chosen for the airspeed measurement based on the AE and the measurement uncorrelatedness assumption.
For the small UAV we are dealing with here, the sampling time Ts = 0.08s is found to be effective for good estimation and detection performance.
5.2.3 MDE and detector operating characteristic curves
A Detector Operating Characteristic (DOC) curve (Sturza, 1988) is a graphical plot that illustrates the power of a discrimination threshold given various fault modes. It is used to understand the trade-off between the false alarm and missed detection rate, and the effectiveness of the designed threshold.
The DOC curves are obtained by plotting various PMD and PFA at different fault sizes. The DOC is the similar to the Receiver Operating Characteristic (ROC) curve, except the y-axis of ROC is the probability of detection PD (PD +PMD = 1). Mathematically, the DOC curve is calculated using Equation (27):
27
where T is the designed threshold and determined by PFA. The functions 
and 
are central and non-central χ2 Cumulative Distribution Functions (CDF), respectively, with Df degrees-of-freedom. The non-centrality parameter λ represents the sum of the non-zero means. In this case, it is sum of the biases vectors 
over the sliding window q, as shown below:
28
The expression for λ here is greatly simplified due to the sequential measurement update procedure. In the next subsection, we will generalize λ expression for dealing with various inhomogeneous measurements under the standard KF measurement update setting.
We use measurements over a short period for λ instead of a single measurement. Not only does this satisfy the Df requirement as mentioned earlier, but also a single faulty measurement may be ineffective. For example, if we wait until the fault grows to 7.5 m/s at 3 seconds, and use this measurement to form the test statistic (e.g., the innovation-based χ2 test statistic only requires one measurement at minimum), then it would be too late in issuing the alarm. Hence we accumulate measurements over a short period and use them to form the detection test statistic.
For a fixed size sliding window q, λ can be calculated using a sliding window of the normalized 
from the past time step k−q + 1 to the current time step k. We define two parameters 
and MDE in Equations (29) and (30):
29
30
where 
is the square root of the sum of normalized fault vectors and MDE the represents the magnitude of the largest fault in the sliding window. The definition of MDE is compatible with the one mentioned in Equation (26) because it represents the maximum tolerable fault size in the sliding window.
Figure 8 shows 
using a window size q = 19 and a constant 
. Each stem represents the sum of the past 19 normalized faulty measurements at sampled time k, where each faulty measurement is simulated based on the linear ramp fault profile shown in the bottom subfigure of Figure 7. The 
represents the sum of the past normalized faulty measurements.
over 3 secs period using the linear ramp fault profile in Figure 7
Figure 9 shows the DOC contour plot for a single degree-of-freedom: Df = mq − n = 1 × 19 − 18 = 1. The plot is created by sweeping different λ over a wide range of PFA and PMD. The contour color is represented by the largest minimum detectable airspeed error 
in the sliding window.
DOC curves where the color represents the associated MDE
If the red dot represents a design choice (PFA, PMD) = (10−5, 10−4), the minimum detectable error requirement MDE is 5.0 m/s, which satisfies the inequality in Equation (26). In other words, the fault detection algorithm can detect the airspeed fault 
after it reaches 5.0 m/s, which corresponds to a detection time τ = 2 s starting from the beginning of the fault profile shown in the bottom sub-figure in Figure 7. Both of those numbers satisfy the given constraint as follows:
31
Of course, this is an ideal situation given the χ2 test statistic is assumed to come from a perfect zero mean, unit variance white sequences. Nevertheless, the analysis provides a systematic way of assessing the fault detection capability for a realistic ramp airspeed fault profile.
Table 2 shows both 
and MDE values over range of PFA and PMD. The numbers in the square bracket are the corresponding MDE values. It is seen that the MDE decreases as PFA and PMD increase and vice versa. This trend is correct and intuitive because a more stringent integrity and continuity requirement would enforce the detection function to catch a fault reliably at a larger MDE since a larger airspeed fault would trigger the detection function to cross the threshold more easily.
The sliding window 
(m/s) and its associated MDE (m/s) in the square bracket for Df = 1
The final sliding window size was chosen to be 19 for both detectors. The sensitivity analysis is done for different sliding window sizes, and it is observed that increasing q from the minimum window size (i.e., 19) does not significantly change the DOC curves. Also, too big of a window size would make the test statistic function sluggish because it tends not to respond slower than the latest change in the measurement.
5.2.4 Trade-off analysis procedure
The above subsections complete the determination of the MDE design based on the integrity requirements (Ireq, Creq) and physical air data system requirement (AE, τmax). Our MDE design can also be used for different pitot tube failure modes (e.g,. stuck or oscillatory fault). Note that analyzing a complete set of pitot tube failure modes would make the MDE calculation statistically more reliable, which is the subject of future studies. In what follows, we summarize the necessary IM performance trade-off design procedures:
Determine a suitable sampling rate for the measurement
Determine a reasonable AE based on the realistic fault mode profile (e.g., constant, ramp, or oscillatory)
Determine a feasible set of requirements Ireq and Creq that satisfies the constraint in Equation (26) using DOC curves
If Ireq and Creq are satisfied, then record the MDE and τ for assessing the fault detection performance
If Ireq and Creq are not satisfied, return step 3. If the pair (Ireq, Creq) is given as a hard requirement, then AE needs to be relaxed (return step 1)
5.2.5 Important trade-off factors
The design procedure is also presented in a flowchart shown in Figure 10. It can be seen that when MDE and τ requirements are not satisfied, either performance requirements Ireq and Creq are needed to be redefined, or AE and τmax are needed to be re-adjusted. We also list the essential trade-off factors and their relationships that need to be considered for the IM performance trade-off design:
Trade-off design procedure
Sampling time versus auto-correlation
Effectiveness of test statistics versus sliding window size
System performance requirements (Ireq, Creq) versus maximum allowable operating conditions (AE, τmax)
Sliding window size versus maximum allowable operating conditions (AE, τmax)
These variables are closely interconnected, and ultimately we need to find a set of (AE,τmax) that is suitable for the given requirements (Ireq, Creq). The choice of (AE,τmax) would largely depend on the size of the sliding window, effectiveness of the chosen detectors (or test statistics), and sampling rate.
The fundamental limitation of change detection is that the design is a compromise between detecting true changes and avoiding false alarms. A poor design gives either a slow filter (no alarm from the test statistic) or a fast filter (many false alarms), and that is the worst can happen.
5.3 MDE analysis for KF residual-based test statistic
The quality of a test statistic is often determined by its MDE (Grosch et al., 2017; Sturza, 1988). As demonstrated in Section 5.2.3, given a set of (PFA,PMD), the MDE (e.g., additive bias) can be found. The MDE informs us of the limits of the fault detection capability. In this section, we show the general expression of MDE for the residual-based test statistic. First, we can express rk−q:k under the faulty H1 hypothesis:
32
where F contains at most q steps since the start of the fault as defined in Equation (A6). F is the fault vector used to represent faults for different types of measurements. If we use both GNSS and airspeed measurement to update KF at the same time, then F can contain at most mq non-zero values. In this case, since the only faulty measurement is assumed to come from the pitot tube, we can further extract the faulty vector and maintain the correct matrix size by applying a column vector E7:
33
where 
. The unit vector e7 is used because the airspeed measurement is on the seventh row shown in Equation (4). The size of I is now 7q × 7q. Then, the expectation of the residual rk−q:k is calculated as follows since 𝔼[W] and 𝔼[V] are zero respectively:
34
Under the faulted hypothesis H1, the residual-based test statistic follows a non-central χ2 distribution with Df equal to mq − n:
35
The non-centrality parameter λ can be computed as:
36
In general, fault can come from any measurement or a combination of different measurements, and the fault characteristic depends on the particular type of measurement used. For an inhomogeneous fault vector F, λ can be expressed as 
without any loss of generality. This might create a challenge if isolation is required. Hence, we continue with the expression shown in Equation (36). The parameter λ can be determined by solving Equation (27) for given a set of PFA, PMD, and Df.
Assuming we have an exact q-step pitot tube faults in the vector F, we can project λ to a q-step detectable error for the pitot tube by reformulating Equation (36):
37
This formulation assumes the fault happens in every time step in the sliding window. In reality, the fault can occur at any time, but it can be only accounted up to q steps. Therefore, the following definition of MDE of 
, denoted as 
, is conservative if q is large, but less conservative when compared to the definition of MDE in Equation (30). The relationship between 
and MDE is shown below:
38
The 
can be interpreted as an average fault over the sliding window and is smaller than the MDE defined in Section 5.2.3. The difference between 
and MDE depends on the sliding window size and the fault profile. If a small sliding window size and a slow ramp fault profile are used, the MDE in the previous subsection is not a bad choice for the fault detection design.
6 PROTECTION LEVEL CALCULATION
In the previous section, we introduced the test statistics and MDE design and analysis for the air data fault detection system. In this section, we derive a new protection level for α and β. We first introduce the definition of the alert limit and protection level in the context of synthetic air data, then we give the formal definition of the protection levels.
6.1 Protection level and alert limit
We define the alert limit of angle-of-attack and sideslip angles needed for the protection level calculation. Alert Limit, denoted as AL, is usually defined as the maximum error in a state estimate that can be tolerated before a system is considered hazardous. Protection Level, denoted as PL, is defined as the guaranteed upper bound of the estimation error uncertainty σ. In theory, we want the probability of the state estimate error ɛ being greater than AL to be extremely low to assure integrity. Practically, we can also formulate the integrity requirement by using PL as follows4:
39
As a consequence, if we have PL < AL, the integrity requirement will be met. The error uncertainty σ is usually inflated by a factor K. This inflated error uncertainty Kσ is the protection level PL.
For a given integrity risk requirement PMD, the PL are calculated. Protection Level PL is a function of the sensing system, and AL is a function of the operation. In other words, PL and AL are independent from each other.
In the case of pitot tube failures being considered here, we are interested in detecting WB faults before they result in the air data estimate (i.e., α and β) being outside of the safe-flight envelope, thereby leading to a control system (or a pilot in the case of a manned aircraft) to execute unnecessary but potentially hazardous maneuvers.
For the purpose of simplicity, we will assume the safe-flight envelope is a rectangle where the upper edge of the safe-flight envelope is defined by the UAVs maximum angle of attack αmax (which is the stall angle of attack αstall). The lower edge is defined by a minimum angle of attack αmin, which is due to some aircraft structural considerations. The left and right edges are defined by βmin and βmax, which are derived either from aerodynamic control or structural strength limits. Therefore, we define the alert limit 
and 
as the absolute nominal safe operating region here since the true reference α and β are typically not available on UAVs. For the UAV considered here, the lower and upper bound of 
and 
is given in the following:
40
These boundaries of the safe-flight envelope form the absolute alert limits for the fault detection algorithm and are generally not the same for different UAVs.
Note that the protection level bound αPL represents a deviation from the true state and the alert limit αAL represents the relative error tolerance in the GNSS applications. However, the absolute alert limit 
in this case is still valid. To see why this is the case, consider the following mathematical equivalence:
41
where αAL represents the alert limit in the relative sense, and the absolute alert limit is defined as 
. Similarly, the absolute protection level is defined as 
. Hence, the new definition of 
does not conflict with the typical definition of αAL. In general, Ireq is a set of discrete probability values representing various fault modes, and the probability of missed detection 
for α, as an example, would be part of the integrity budget Ireq. Nevertheless, we will use the same value for the protection level calculation since we are only dealing with one fault mode, that is, 
.
Figure 11 graphically depicts both protection level 
and alert limit 
in relation to α. Ideally, 
should rarely go over 
due to the small integrity risk requirement. When 
does exceed 
, we can safely conclude it is highly likely to have been the result of a faulty pitot tube and the integrity of α is lost.
Depiction of protection level and alert limit of α
6.2 Protection level for synthetic air data
The fault detection algorithm needs to guarantee (at a certain level of confidence) that estimation error in α and β do not exceed their alert limits. Various slope-based PL calculations have been used to protect horizontal and vertical state errors against single (Walter & Enge, 1995; Brown & Chin, 1998; Milner & Ochieng, 2011) or multiple (Pervan et al., 1998; Angus, 2006; Blanch et al., 2009; Jiang & Wang, 2014) GNSS faults. However, these methods are developed for solving the redundant measurement problem.
Furthermore, the PL is usually calculated to protect states in an EKF. We develop a PL method that can protect the states derived from the EKF states. In particular, we calculate the PL for the synthetic angle-of-attack α and sideslip β estimates.
Under the fault-free hypothesis H0, the PL is calculated by inflating state errors by a factor of K. In particular, the PL of α and β are calculated as follows:
42
43
where Kα,0 and Kβ,0 under H0 are calculated as follows:
44
where Q is the tail distribution function of the standard normal cdf. The matrix Pk is the state covariance from the EKF and Aαβ is the flow angle propagation transformation matrix specified in Sun et al. (2019b). The unit vectors e1 = [1, 0]T and e2 = [0, 1]T are used to extract the diagonal terms of 
.
In the presence of fault under the hypothesis H1, PL needs to be increased to account for the faulty airspeed component 
. The formulation is done as follows according to (Brown & Chin, 1998; Angus, 2006):
45
46
where slopeα and slopeβ represent the ratio between the α or β state error and the standard deviation of the test statistic Dr,k. The notion of the maximum slope (Brown & Chin, 1998) is not applicable here since there is no redundant airspeed at each time step to calculate a set of slopes. Specifically, since the accumulated sequential residual or innovation vectors come from the single pitot tube, so the source of the fault is always known and singular. In other words, the multiple hypothesis tests of determining which GNSS measurement is faulty by calculating the maximum slope does not apply here. The slope value in this study depends on the test statistic window size and severity of the fault profile.
The inflation factors Kα,1 and Kβ,1 under H1 can be calculated as follows:
47
where 
represents the failure rate of the pitot tube due to water blockage. Generally, the failure rate of the pitot tube due to a particular fault is calculated based on rigorous sensor characterization testing. Since the failure rate testing is not done for this study, we assume K(·),1 = K(·),0 = Q−1(PMD/2) for both α and β. This assumption is valid but conservative as K(·),1 is generally greater than K(·),0 in general. The parameter λU is the upper confidence bound for the maximum q-step airspeed faults, and it is computed as follows:
48
This upper bound λU can make the protection levels overly conservative if the size of the sliding window q is large. Hence, the tightness of the protection level is another factor for choosing an appropriate q.
It is worth noting that many small UAVs do not have angle vane sensors to provide a set of α and β. The protection levels of α and β are particularly useful when the angle vane sensor is not available. We can monitor α and β based on the synthetic air data estimates and protect the vehicle from exceeding α and β flight envelopes due to the pitot tube failure.
7 FLIGHT DATA TESTING
The fault detection algorithm is tested using a flight data set recorded by a UAV. The UAV is the Sentera Phoenix (PHX) shown in Figure 12. It has an Eagle Tree pitot tube (Eagle Tree Systems, 2020) attached to the right wing to measure airspeed. PHX utilizes a similar version of Pixracer autopilot (Pixhawk FMU-V4) for control and navigation.
PHX by Sentera LLC
For this particular flight operation, the UAV crashed 45 seconds into the flight due to a water-plugged faulty pitot tube. This flight data is suitable for the fault detection algorithm analysis since it contains a known WB pitot tube fault signature.
Also, since the UAV has only one pitot tube, only one SADS is employed to show the flight results. Ideally if this UAV carried two pitot tubes and two independent SADS, then it would follow the decision logic in Table 1 to raise alarm and switch to SADS-2 when SADS-1 detected and isolated the fault airspeed measurement.
Figure 13(a) and 13(b) show the trajectory, altitude, and airspeed over the short flight period before the crash. The UAV takes off around 570 s and circles up to about 70 m, then flies straight northwest direction to the edge of the crop field. It is seen in Figure 13b that the airspeed of UAV experiences a sharp drop right after 618 s. The airspeed measurement eventually drops to a negative value. The last recorded altitude is around 35 m by the onboard autopilot.
Trajectory and airspeed information over the time of flight
Figure 14 shows the attitude estimates from the SADS filter. Notice that the Euler angle estimates start to change abruptly at the end of the flight due to the faulty pitot tube. For example, the pitch angle θ and roll angle ϕ increased dramatically around 618 s. The high value θ indicates the UAV might have been pitching up and stalling.
SADS attitude estimates over time
Figure 15 shows three different test statistics and their corresponding thresholds from the actual takeoff timestamp of 570 s to the end of the flight.
Detection variables over time
The top sub-figure shows the sliding window innovation-based test statistic Dγ over time where the window size q is 19. The bottom sub-figure shows a residual-based test statistic Dr and a Geometric Moving Average (GMA) residual-based test statistic 
, where the window size q is also 19 for both. The 
is defined as the same as Dr except the residual vector 
is written as follows:
49
In Equation (49), it is seen that the 
uses an exponential forgetting factor μi to weigh on the past measurements less than the recent. The innovation-based threshold Tγ and residual-based Tr are different from each other due to the difference in Df, even though the same size of the sliding window is used.
It is seen that the innovation-based test statistic exceeds its threshold three times. It appears that the first occurrence at 593 s is most likely a false alarm given the fact we know the pitot tube WB fault occurs at the end of the flight. Furthermore, by examining the attitude in Figure 14, the large ϕ and sudden change in θ suggests a momentary accelerated stall.
It is well known that the innovation-based test statistic is sensitive to the highly dynamic motion. The UAV experiences a large roll change at 593 s, which may cause Dγ to think there might be a fault in the system. The second occurrence might come from a sudden increase in airspeed at 610 s. Though we do not know if the fault has occurred yet by visually examining the airspeed plot, Dγ indicates there might be a fault.
Notice Dγ goes below the threshold at 614 s before rising to cross the threshold again at 618 s for the third time. This unstable behavior is not a good characteristic for a detector because it will be issuing many false alarms and failing to provide a crisp decision of a fault’s occurrence.
It is also well known that the innovation-based test statistic is sensitive to sensor noise scaling (Gustafsson, 2000), and the innovation covariance S inside of Dγ is sensitive to noise scaling. The relationship is illustrated by Equation (50): S is affected by noise covariance matrices Rw and R, and the initial state covariance P0. Small η can make Dγ sensitive to noise even though the state estimates are not affected.
50
Though not shown here, by scaling the process noise 
of the wind vector without adversely changing the performance of the estimator, the first two crossing occurrences of Dγ can be suppressed. However, hand-tuning is not recommended because it might adversely affect other estimates.
On the other hand, both residual-based test statistics Dr and 
exceed their corresponding thresholds only once at the end of the flight. The threshold value for Dr and 
are the same, hence only one Tr is plotted. Though small increases at 593 s and 612 s can be seen in Dr, Dr did not cross its threshold.
The residual-based test statistic is less sensitive to both highly dynamic motion and noise scaling due to the larger weighting factor Σ. Furthermore, the GMA test statistic 
appears to be even less sensitive. 
is discounting many early measurements in the sliding window. Also, it rises faster than Dr when crossing the threshold, indicating a shorter fault detection time. When comparing detection time alone, the innovation does appear to have the fastest detection at the third crossing.
The GMA technique used here illustrates that sometimes different techniques can be used to improve the baseline detector function. For example, we may combine the CUSUM and residual-based test statistics to improve detectability. Ultimately, it is a trade-off between the detection time and false alarm when choosing a detection test statistic. In this case, it appears the GMA residual-based test statistic 
is the best for the faulty pitot tube detection since it crosses the threshold at the correct incident and provides a good detection time.
In order to see how observability can actually play an important role in fault detection, we also monitor the observability Gramian of the system over time. Different metrics can measure the degree of discrete observability Gramian, such as the determinant, trace, or the condition number (Summers et al., 2016; Avant & Morgansen, 2019). In work here, we utilize the condition number κ to quantify the degree of observability. That is:
51
where 𝒢d,k−q:k the discrete observability Gramian using the information from the time step k-q to the current k. In general, if the condition number of the observability Gramian is large, it means the states are not well observed. The observability Gramian here only requires a finite horizon instead of the infinite horizon. This calculation is done to monitor the recent motion of the dynamics of the vehicle. It also reduces the computational burden for the computer processor.
Figure 16 shows the normalized condition number of observability Gramian over time. The poor condition number is expected before takeoff since the observability Gramian rank is deficient. The condition number quickly reaches near zero after 570 s (takeoff). However, it is seen that the condition number got a bit worse from 600 s to 607 s. This change is perhaps due to the straight level flight (no heading change), or some water started to enter the pitot tube before it became evident.
Normalized condition number of observability Gramian
The second claim is deduced based on examining all the EKF state estimates at 600 s. Since there is no abnormal phenomenon from the position, velocity, and attitude estimates around 600 s, the IMU and GNSS measurements are assumed to be nominal. Hence, the culprit is either the poor observability due to the trajectory or the pitot tube.
A poor condition number can potentially trigger the fault detection algorithm to raise the alarm even though there is no real fault. If there is a fault that is about to happen, then monitoring the observability Gramian can potentially warn the system before the fault occurs. If a poor condition number is a result of the straight-level flight, then weak observability can mask the fault detectability if a fault occurs during this time. Therefore, a close examination should be carried out when using observability to detect the actual fault.
Figures 17(a) and 17(b) show the α and β estimates, 2σ uncertainty bounds, protection levels αPL and βPL, and alert limits αAL and βAL. The protection level bound is for the single SAD. The protection level αPL goes outside the alert limit momentarily at 590 s, 592 s, and 595 s, and exceeds both lower and upper alert limit at the end of the flight. The s-turn causes the first three crossings before heading to northwest direction. The last one is caused by the faulty pitot tube. It can be seen that the α increases drastically at the end of the flight, which leads to the stall and eventual crashing. Sideslip β also experiences a similar change during the s-turn. The sideslip estimate β changes from positive to negative, then positive again around 590 s, which is intuitively correct based on the s-turn and the northeast tail wind direction.
Flow angle α and sideslip β estimate, uncertainty bound, protection level, and alert limit
The protection level βPL exceeds the alert limit a couple of times and eventually goes outside the alert limit at the end of the flight as expected. The protection level of α and β can be used to check the validity of ADS integrity, which is much more useful than just looking at the confidence uncertainty bound represented by 2σ. For example, we can see that the α and its 2σα seem to be acceptable from 593 s to 596 s, but the integrity of the α estimate is actually lost during this time based on the protection level. This loss of integrity can be used as another flag for the validity of the α estimate.
8 CONCLUSION AND FUTURE OUTLOOK
A dual pitot tube ADS is designed with fault detection and isolation capability for small UAS. The purpose of this algorithm is to provide a reliable ADS and recovery strategy for safe drone operations in case of the pitot tube water-blockage fault under rainy and foggy conditions. The fault detection algorithm is designed to detect faults based on the given integrity requirements and known water-blockage fault profile.
Systematic procedures and various factors are laid out to show how to design the fault detection algorithm. Though the stringent performance requirement, such as 10−7, may not be realizable to certify ADS in small UAS due to the low-cost sensors onboard, this IM approach allows engineers to assess the performance of the FDI algorithm from the requirement point of view. The high-performance requirements can potentially be achieved for the ADS in small UAS when highly accurate sensors (e.g., good GPS and IMU), good control design (reject external disturbances), and sensible path planning (observability-aware trajectory design) are employed.
HOW TO CITE THIS ARTICLE
Sun, K., & Gebre-Egziabher, D. (2021). Air data fault detection and isolation for small UAS using integrity monitoring framework. NAVIGATION, 68, 577–600. https://doi.org/10.1002/navi.440
ACKNOWLEDGMENTS
The authors gratefully acknowledge the Minnesota Invasive Terrestrial Plants and Pests Center (MITPPC) through the Minnesota Environment and Natural Resources Trust Fund for financial support to conduct research associated with increasing the reliability of small UAV technology used for surveying applications. The authors also gratefully acknowledge Todd Colten and Sentera LLC for donating the flight data for the air data fault detection research.
APPENDIX A: BATCH LINEAR SYSTEM REALIZATION
In this appendix, we present one form of the batch linear system realization to faciliate the connection between the LTV observability matrix and the residual-based test statistic shown in Section 5. We start with a model of the discrete linear dynamic system with an unknown additive fault:
A1
where fk is the deterministic additive fault vector, and xk, zk, υk and υk are defined in Section 2.2. The vector wk is the process noise vector and is assumed to follow N(O, Rw). The matrices Φk and Γk are the state transition matrix and discrete noise coefficient matrix, respectively.
Only additive faults are chosen for this study to analyze the feasibility of air data fault detection capability. Nonlinear fault with more sophisticated models should be considered in the future.
Since there are not enough redundant airspeed measurements at each time step k, we obtain the following batch realization by stacking all the measurement vectors from the past time k−q to the current time step k in Equation (A2):
A2
The matrices Zk−q:k, Wk−q:k, Vk−q:k, Fk−q:k, 
and Qw,k−q:k are defined in the following:
A3
A4
A5
A6
A7
A8
The matrix 
is the discrete LTV observability matrix over a sliding window. This batch realization [extension of the linear time-invariant system in (Isermann, 2005)] formulation connects the stacked measurements Zk−q:k to the past state vector xk−q and its associated observability matrix 
nicely.
Batch realization is not unique; Joerger and Pervan (2013) define a different formulation where the past measurements are a function of all the past states. However, the measurement model in Joerger and Pervan (2013) is not explicitly expressed as a function of the observability matrix.
APPENDIX B: PROOF OF MATRIX EQUALITY
In this appendix, we prove the following equation is true:
B1
where 
is the observability matrix shown in Equation (A7), 
is left pseudoinverse of 
shown in Equation (21), and Σ is the weighting matrix shown in Equation (18). Note that Σ is symmetric by construction, that is, Σ = ΣT
Proof. First, we show 
is symmetric:
B2
where the 7th equality is achieved because 
is also symmetric.
Now we have the following fact: (B3)
B3
With the equality above, we can express Equation (B1) as follows:
B4
Now we proceed to prove the following: 
.
This can be shown to be true by moving the left-hand side of equation to the right as follows:
B5
We also realize the following equation is true:
B6
Note that 
is justified because 
is right-invertible. However, the right inverse (not Moore-Penrose) of 
is not unique.
By substituting Equation (B6) into Equation (B5), we show that 
is indeed true. Finally, we arrive the matrix property that we want to prove:
B7
Footnotes
Funding information
Minnesota Invasive Terrestrial Plants and Pests Center (MITPPC)
1 The temperature and altitude effect are not considered in this study
2 The second airspeed data shown here is digitally extracted from the paper for illustration
3 Since a sliding window of measurement is used, the test statistics always has enough nominal measurements for χ2 test before the fault occurs. If the faulty measurement enters the sliding window too slowly, the test statistics can be ineffective
4 This expression is simplified from the following: P(|ɛ| > PL, D < T|Hi)P(Hi) = P(|ɛ|> PL|Hi)P(D < T|Hi)P(Hi) ≤ P(|ɛ| > PL|Hi) P(Hi) ≤ PMD, where D is a detector function or test statistic, and T is a designed threshold
- Received October 2, 2020.
- Revision received May 10, 2021.
- Accepted July 6, 2021.
- © 2021 Institute of Navigation
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
- Article
- Abstract
- 1 INTRODUCTION
- 2 DUAL PITOT TUBE AIR DATA SYSTEM DESIGN
- 3 AIR DATA SYSTEM REQUIREMENT
- 4 WATER BLOCKAGE FAILURE MODE
- 5 FAULT DETECTION DESIGN AND ANALYSIS
- 6 PROTECTION LEVEL CALCULATION
- 7 FLIGHT DATA TESTING
- 8 CONCLUSION AND FUTURE OUTLOOK
- HOW TO CITE THIS ARTICLE
- ACKNOWLEDGMENTS
- APPENDIX A: BATCH LINEAR SYSTEM REALIZATION
- APPENDIX B: PROOF OF MATRIX EQUALITY
- Footnotes
- REFERENCES
- Figures & Data
- Supplemental
- References
- Info & Metrics
Related Articles
Cited By...
- No citing articles found.