## 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 *V*_{a}, 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.

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.

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.

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.

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, *δ***v**^{n} = [*δV*_{N} *δV*_{E} *δV*_{D}]^{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 *δ***b**_{a} = [*δb*_{ax} *δb*_{ay} *δb*_{az}]^{T} and *δ***b**_{g} = [*δb*_{gx} *δb*_{gy} *δb*_{gz}]^{T} are the accelerometer and rate gyro triad output bias error vectors, respectively. Finally, *δ***W**^{n} = [*δW*_{N} *δW*_{E} *δW*_{D}]^{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 **v**^{n} estimates from the GNSS receiver, along with the true airspeed *V*_{a} 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 **R**_{w} 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 *δ***z**_{k} = **H**_{k}*δ***x**_{k} + **υ**_{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 **H**_{k} ∈ ℝ^{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

**W**^{n}must be*quasi-static*. The term*quasi-static*means that the variations in**W**^{n}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 **W**^{n} 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 *I*_{req} and continuity risk *C*_{req}. 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 *P*_{MD} represents the probability of not detecting a pitot tube failure given the pitot tube has indeed failed. Similarly, the probability of false alarm *P*_{FA} 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 *H*_{0}. If the pitot tube is not working correctly, this hypothesis is denoted as *H*_{1}. 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 *I*_{req} and *C*_{req} 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.

The airspeed *V*_{a} is typically calculated based on Bernoulli’s principle as follows^{1}:

13

where *P*_{t} is the stagnation or total pressure, *P*_{s} 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.

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 here^{2} 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*/*s*^{2} and 3 *m*/*s*^{2}, 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 **v**^{n} and wind vector **W**^{n} 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.

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 *T*_{s}, the sliding window size *q* of the fault test statistic, the integrity risk *P*_{MD}, and continuity risk requirements *P*_{FA}. 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 *H*_{0} as shown below:

15

where *D*_{γ, k} is test statistic at the time step *k*, *γ*_{k} is the innovation vector calculated using , the matrix *S*_{k} 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 **r**_{k−q:k} is stacked residual vector from the past time step *k*−*q* to the current time step *k*, denoted as . Each residual **r**_{k} is computed from the difference between the measurements and a posteriori estimate . The matrix **Σ** is the covariance matrix of the weighted least residual vector **r**_{k−q:k} is shown in Equation (18):

18

where ⊗ is the Kronecker tensor product and **R**_{w} 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 *x*_{k−q} using weighted least squares. The state vector *x*_{k−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 **r**_{k−q:k} can be subsequently expressed in the following:

22

Under the fault-free *H*_{0} hypothesis, we can also write **r**_{k−q:k} as follows:

23

where **r**_{k−q:k} follows a normal distribution *N*(**O**, **Σ**), and **W**_{k−q:k} and **V**_{k−q:k} are defined in Equation (A4) and (A5) respectively.

Using Equation (22), we can also write *D*_{r,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 *D*_{r,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 *H*_{0} 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 *V*_{a,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 *V*_{a,stall}.

Since the UAV’s operating speed *V*_{a} 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 *V*_{a,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*/*s*^{2}. 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 *V*_{a,stall} translates to the maximum detection time *τ*_{max} = 3*s* 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 *P*_{FA} and *P*_{MD}. 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*/*s*^{2}. The fault detector should still catch any rate that is higher 2.5 *m*/*s*^{2} 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 *T*_{s} = 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}.

On the other hand, the sampling time *T*_{s} = 0.08*s* is sufficient even if we only use the measurements after the fault occurs to form the test statistic^{3}. 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 *T*_{s} = 0.08*s* 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 *P*_{MD} and *P*_{FA} 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 *P*_{D} (*P*_{D} +*P*_{MD} = 1). Mathematically, the DOC curve is calculated using Equation (27):

27

where *T* is the designed threshold and determined by *P*_{FA}. 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.

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 *P*_{FA} and *P*_{MD}. The contour color is represented by the largest minimum detectable airspeed error in the sliding window.

If the red dot represents a design choice (*P*_{FA}, *P*_{MD}) = (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 *P*_{FA} and *P*_{MD}. The numbers in the square bracket are the corresponding MDE values. It is seen that the MDE decreases as *P*_{FA} and *P*_{MD} 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 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 (*I*_{req}, *C*_{req}) 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

*I*_{req}and*C*_{req}that satisfies the constraint in Equation (26) using DOC curvesIf

*I*_{req}and*C*_{req}are satisfied, then record the MDE and*τ*for assessing the fault detection performanceIf

*I*_{req}and*C*_{req}are not satisfied, return step 3. If the pair (*I*_{req},*C*_{req}) 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 *I*_{req} and *C*_{req} 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:

Sampling time versus auto-correlation

Effectiveness of test statistics versus sliding window size

System performance requirements (

*I*_{req},*C*_{req}) 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 (*I*_{req}, *C*_{req}). 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 (*P*_{FA},*P*_{MD}), 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 **r**_{k−q:k} under the faulty *H*_{1} 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 **E**_{7}:

33

where . The unit vector **e**_{7} is used because the airspeed measurement is on the seventh row shown in Equation (4). The size of **I** is now 7*q* × 7*q*. Then, the expectation of the residual **r**_{k−q:k} is calculated as follows since 𝔼[**W**] and 𝔼[**V**] are zero respectively:

34

Under the faulted hypothesis *H*_{1}, 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 *P*_{FA}, *P*_{MD}, 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 follows^{4}:

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 *P*_{MD}, 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, *I*_{req} 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 *I*_{req}. 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.

### 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 *H*_{0}, 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 *H*_{0} are calculated as follows:

44

where *Q* is the tail distribution function of the standard normal cdf. The matrix *P*_{k} 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 **e**_{1} = [1, 0]^{T} and **e**_{2} = [0, 1]^{T} are used to extract the diagonal terms of .

In the presence of fault under the hypothesis *H*_{1}, 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 *D*_{r,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 *H*_{1} 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}(*P*_{MD}/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.

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.

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.

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.

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 *D*_{r} 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 *D*_{r} 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 *T*_{r} 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 **R**_{w} and **R**, and the initial state covariance **P**_{0}. 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 *D*_{r} and exceed their corresponding thresholds only once at the end of the flight. The threshold value for *D*_{r} and are the same, hence only one *T*_{r} is plotted. Though small increases at 593 s and 612 s can be seen in *D*_{r}, *D*_{r} 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 *D*_{r} 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.

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.

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 **f**_{k} is the deterministic additive fault vector, and **x**_{k}, **z**_{k}, **υ**_{k} and **υ**_{k} are defined in Section 2.2. The vector **w**_{k} is the process noise vector and is assumed to follow *N*(**O**, **R**_{w}). 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 **Z**_{k−q:k}, **W**_{k−q:k}, **V**_{k−q:k}, **F**_{k−q:k}, and **Q**_{w,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 **Z**_{k−q:k} to the past state vector *x*_{k−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 7^{th} 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 ineffective4 This expression is simplified from the following:

*P*(|*ɛ*| >*PL*,*D*<*T*|*H*_{i})*P*(*H*_{i}) =*P*(|*ɛ*|>*PL*|*H*_{i})*P*(*D*<*T*|*H*_{i})*P*(*H*_{i}) ≤*P*(|*ɛ*| >*PL*|*H*_{i})*P*(*H*_{i}) ≤*P*_{MD}, 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.