Drone Ground Impact Footprints with Importance Sampling: Estimation and Sensitivity Analysis ^†

Abstract

This paper addresses the estimation of accurate extreme ground impact footprints and probabilistic maps due to a total loss of control of fixed-wing unmanned aerial vehicles after a main engine failure. In this paper, we focus on the ground impact footprints that contains 95%, 99% and 99.9% of the drone impacts. These regions are defined here with density minimum volume sets and may be estimated by Monte Carlo methods. As Monte Carlo approaches lead to an underestimation of extreme ground impact footprints, we consider in this article multiple importance sampling to evaluate them. Then, we perform a reliability oriented sensitivity analysis, to estimate the most influential uncertain parameters on the ground impact position. We show the results of these estimations on a realistic drone flight scenario.

1. Introduction

Assessing the risks and feasibility of unmanned aerial vehicle (UAV) operations for outdoor inspection or monitoring missions has become a major challenge for regulatory authorities and drone operators. This evaluation relies on risk analysis methods that can be helpful in the process of flight authorization, but also in the design and the preparation of the mission. Two main types of methods are classically used. The first one relies on the qualitative evaluation of risks by applying some predefined methodologies or guides [1]. This is, for example, the case of classical methods such as failure modes and effects analysis (FMEA) or, more recently and more specifically developed for UAV operations, SORA (specific operation risk assessment) [2]. The second type of methods relies on the quantitative evaluation of risks, based on the use of models developed to represent the UAV behavior, its environment, etc. This is the case of model-based probabilistic risk assessment (PRA) approaches that have recently gained a huge interest for UAVs, see e.g., [1,2,3,4,5]. With these approaches, the accuracy of models used for risk probabilities’ computations is of paramount importance. Indeed, being too conservative may prevent or restrict some operational uses of UAVs, while not being conservative enough may lead to uses with uncontrolled risks. A fundamental keystone in these methods, when considering ground risk evaluation, is the computation of probabilities of impact of the UAV at ground level. Accurate models should be developed to be able to compute representative predictions of impact points’ locations and probabilities. Works from the literature have focused on computing impact point locations, enabling to obtain estimates of impact footprints on the ground level. In [6], impact footprints are computed by reachability analysis, considering a gliding descent model for a fixed-wing UAV, composed of a turning and a straight line phase. Different modes of failure (engine, engine+rudder+ailerons) are considered, as well as the effect of the altitude of the vehicle at the failure instant. This type of impact footprints has also been obtained in [3], considering a 6 degrees-of-freedom dynamic model of a fixed-wing aircraft. The effect of wind on impact footprints has been investigated in [7]. Computation of impact locations and footprints has also been performed in [8], for both fixed-wing and multi-rotor UAVs considering different modes of failure. Some level sets are computed to provide some insights on the distribution of the impact points inside the footprint.

The generation of probability maps has been investigated in other works to provide more information on the impact distributions on the ground that could be useful and reduce conservatism in risk analysis or decision making. In [9], a ballistic model with drag force is considered to represent the descent of a fixed-wing UAV and generate probability density functions of impact points. Uncertainties on drag, initial speed at the instant of failure and external wind are accounted for. Full flight dynamics of a Cesna 182 aircraft are used in [10] to compute ground impact probability maps by Monte Carlo simulations. Total loss of power is assumed and uncertainties on the initial conditions of the UAV at failure instant as well as on the deflection of unactuated control surfaces are considered. A 6 degrees-of-freedom flight mechanics model is also used in [11] for a fixed wing UAV to estimate ground impact probability maps, taking into account the influence of wind direction and speed. Real flight data have been used to model uncertainties on the turn rate and flight path angle of the vehicles for cruise-like mode at constant altitude and straight line. These uncertainties along with the ones on the actuators deflections at the instant of failure are used in the Monte Carlo process. Influence of initial altitude, speed and wind (speed and direction) are analyzed, and a full data basis has been obtained containing impact probability maps for a sampled set of values for these quantities. This data basis can be useful for risk evaluation along a given UAV flight trajectory, e.g., for mission preparation.

Since Monte Carlo simulations can be time-consuming, more recent works have been dedicated to the development of surrogate models for the generation of ground impact probability maps. K-Nearest neighbors models have been considered in [10] to approximate impact probability distribution. Other techniques such as Krigging have been investigated [7] regarding impact footprints or neural networks for both generation of impact footprints [7] and probability maps [11].

In all these works, assumptions are made on the uncertainties on the variables used as inputs of the computations. Uncertainty representations are mainly based on statistic models (probability distributions) and/or bounds (intervals with no statistical assumptions). Accuracy of the resulting outputs (ground impact locations, footprints, probability maps) strongly relies on the representativeness of these assumptions. Another important aspect in these approaches is the computation method itself and the choice of its hyper-parameters. For example, choice of simulation budgets in Monte Carlo approaches is crucial, as it may strongly influence the probability density estimation and its confidence.

Moreover, Monte Carlo methods with a low number of samples lead to an underestimation of extreme ground impact footprints, which may be of interest to provide more confidence in the risk assessment process for flight preparation and authorization. Knowledge of UAV’s extreme fallout zones can also help defining safety levels at very low thresholds, which can be critical for certain high-risk infrastructures.

This paper therefore addresses the estimation problem of accurate extreme ground impact probability maps and footprints containing 95%, 99% and 99.9% of the impacts. Multiple importance sampling (MIS) is considered to estimate density minimum volume sets associated with these extreme quantiles.

In addition, a study of the sensitivity of hazard parameters is proposed to estimate the most influential uncertain parameters on ground impact positions. This analysis may enable both operators and drone constructors to better understand, design and anticipate fallout zones in the event of a failure. All the results in this paper are obtained for the case of a fixed wing UAV after main engine failure.

This paper is organized as follows. The first section focuses on the simulation approach used to compute the ground impact points coordinates. The characterization of uncertainties that are taken into account is discussed in Section 3. The estimation of ground footprints by Multiple Importance Sampling is then presented in Section 4 and an associated sensitivity analysis is proposed in Section 5.

2. Ground Impact Simulation

In this paper, we focus on impacts on the ground due to a loss of control of the UAV (unmanned aerial vehicle) after a main engine failure. It is assumed that immediately after the failure, the engine thrust becomes zero and the control surfaces remain stuck in their equilibrium positions. The objective of this section is to present the models and approach that are used to compute the impact points at ground by simulation, based on previous studies by the authors in [7].

2.1. UAV Dynamics

The model used to simulate the trajectory of the UAV to the ground is a six degrees of freedom (6DOF) dynamic model, including full flight mechanics, and hence enabling one to incorporate the influence of wind from a dynamical point of view (and not kinematical compared to some approaches developed in the literature [12]). The model considered here is a fixed wing aircraft such as the one presented in [13]. The control input vector $u={\left[\begin{array}{cccc}{\delta }_a& {\delta }_e& {\delta }_r& {\delta }_T\end{array}\right]}^{\top }$ is composed of ailerons, elevators, rudder deflections, and thrust command. The state of the dynamical system to be simulated is defined as $\chi ={\left[\begin{array}{cccc}{X}^{\top }& V^{\top }& {\eta }^{\top }& {\Omega }^{\top }\end{array}\right]}^{\top }$ where $X={\left[\begin{array}{ccc}x& y& z\end{array}\right]}^{\top }$ is the position vector defined in a local NED (north east down) frame, V and $\Omega$ are the translation and angular velocity vectors in the aircraft body-frame, and $\eta ={\left[\begin{array}{ccc}\varphi & \theta & \psi \end{array}\right]}^{\top }$ is the vector of Euler angles (roll-pitch-yaw) describing the attitude of the UAV. The origin of the (inertial) local NED frame is chosen to correspond to z = 0 (ground) and is arbitrarily chosen for x and y-components, as we are only interested in the description of the motion of UAV during its descent to the ground with respect to the vehicle position at the instant of failure (considered to be $x=y=0$).

The rigid-body dynamics of the UAV is described as

$$\left\{\begin{array}{c}\dot{X}=R_{\eta }V\hfill \\ \dot{V}=-\Omega \times V+\frac{1}{m}F\hfill \\ \dot{\eta }=T_{\eta }\Omega \hfill \\ \dot{\Omega }=J^{-1}(-\Omega \times J\Omega +M)\hfill \end{array}\right.$$

(1)

where $R_{\eta }$ is the orientation matrix parametrized in terms of Z-Y-X Euler angles given by

$$R_{\eta }=\left[\begin{array}{ccc}{c}_{\theta }{c}_{\psi }& s_{\varphi }{s}_{\theta }{c}_{\psi }-s_{\psi }{c}_{\psi }& s_{\varphi }{s}_{\psi }+s_{\theta }{c}_{\varphi }{c}_{\psi }\\ c_{\theta }{s}_{\psi }& c_{\varphi }{c}_{\psi }+s_{\theta }{s}_{\varphi }{s}_{\psi }& s_{\theta }{c}_{\varphi }{s}_{\psi }-s_{\varphi }{c}_{\psi }\\ -s_{\theta }& c_{\theta }{s}_{\varphi }& c_{\theta }{c}_{\varphi }\end{array}\right]$$

(2)

and $T_{\eta }$ is the transformation matrix defined by

$$T_{\eta }=\left[\begin{array}{ccc}1& s_{\varphi }{t}_{\theta }& c_{\varphi }{t}_{\theta }\\ 0& c_{\varphi }& -s_{\varphi }\\ 0& s_{\varphi }/c_{\theta }& c_{\varphi }/c_{\theta }\end{array}\right]$$

(3)

whith the notations $c_{\alpha }=cos\left(\alpha \right)$, $s_{\alpha }=sin\left(\alpha \right)$ and $t_{\alpha }=tan\left(\alpha \right)$ for any given angle $\alpha$. The inertia matrix of the UAV is denoted by J and its mass by m. Values used for the UAV parameters are given in Appendix A.

The resulting force F expressed in the aircraft body-frame

$$F=F_{eng}\left({\delta }_T\right)+F_g\left(\eta \right)+F_a(V,V_w,\eta ,\Omega ,{\delta }_a,{\delta }_e,{\delta }_r)$$

(4)

is composed of the thrust, due to the engine ($F_{eng}\left({\delta }_T\right)$) which is zero after engine failure, the gravity force ($F_g\left(\eta \right)$) and the resulting aerodynamic force $F_a$, which depends on the true air speed of the UAV and then on the wind speed vector $V_w$. To compute this force, a full aerodynamic model is used, such as the one described in [13] and is presented in Appendix A. Similarly, the resulting torque expressed in the aircraft body-frame

$$M=M_{eng}\left({\delta }_T\right)+M_a(V,V_w,\eta ,\Omega ,{\delta }_a,{\delta }_e,{\delta }_r)$$

(5)

is composed of the torque component $M_{eng}\left({\delta }_T\right)$, due to the thrust of the main engine (equals to zero after engine failure) and the aerodynamic torque $M_a$, which also depends on the wind speed vector $V_w$. To compute this torque, a full aerodynamic model of the aircraft is also considered (see Appendix A). The dynamic model of the UAV can be summarized with the following state-space representation

$$\dot{\chi }=f(\chi ,u,V_w)$$

(6)

The simulation of the UAV descent trajectory is performed from an initial condition ${\chi }_0$, defined at the engine failure instant $t_0$, to ground impact, that corresponds to instant $t_f$ such that the altitude $h\left(t_f\right)=-z\left(t_f\right)=0$.

During a steady flight (coordinated turn, straight flight, pull-up/pull-over etc.), the accessible space through the initial condition $({\chi }_0,u_0)$ is considerably reduced. Defining the initial condition then consists of zeroing numerically the dynamic part of Equation (6), while simultaneously considering kinematic constraints related to flight mode [14]. In this case, the control vectors and the dynamic part of the state vector are entirely defined by these constraints. This method is called trim algorithm. A simple way to represent a trajectory is to consider the two following parameters:

• the turn rate $R=d\psi /dt$, where $\psi$ is the heading angle
• the flight path angle $\gamma =\dot{z}/V_a$, where $V_a$ is the aerodynamic speed of the aircraft.

Therefore, the trim algorithm can be run to determine the initial condition $({\chi }_0,u_0)$, by assigning values $R_0,{\gamma }_0,V_{a_0}$ and $h_0$ to the turn rate, flight path angle, aerodynamic speed and altitude, which are representative of the UAV flight conditions. A straight cruise flight mode at constant altitude can, for example, be considered by choosing $R_0=0$ and ${\gamma }_0=0$. Note that the wind is not considered in the trim algorithm.

From this initial condition $({\chi }_0,u_0)$, the trajectory of the UAV is simulated until the impact time $t_f$, by considering zero thrust (main engine failure) and taking into account the wind speed $V_w$. The complete simulation process is represented in Figure 1.

Using this approach, one can simulate a single trajectory and compute the coordinates of the ground impact point. An example of the trajectory is provided on Figure 2, corresponding to ${\gamma }_0=0$ deg, $R_0=0.15$ rad/s, $V_{a_0}=30$ m/s, $h_0=150$ m and no wind.

Nevertheless, impact point coordinates are not deterministic in the sense that their computation will suffer from different sources of uncertainties on the parameters of the problem. The next section covers the characterization of these uncertainties, as well as the Monte Carlo approach to handle them.

3. Uncertainties and Monte Carlo approach

The computation of an impact point involves a simulation relying on several parameters. Some of them will be considered as fixed values, such as the initial ground altitude $h_0$ and aerodynamic speed $V_{a_0}$. Note that these quantities are affected by some uncertainties, since the reference altitude and velocities commanded for the UAV are not exactly flown in practice. However, the influence of their respective incertitude levels on the impact point location is negligible. Uncertainties on the parameters of the UAV model (aerodynamic coefficients) are also not taken into account in this paper. A robustness analysis with regard to them should be carried out, especially since the descent phase to ground may lead to aerodynamic behaviors different from the ones that can be identified. This is beyond the scope of this paper and will be considered in future work. Uncertainties taken into account in this article are described in the following subsections.

3.1. Uncertainties on $R_0$ and ${\gamma }_0$

Experimental data have been recorded on flights realized by Altametris, the drone subsidiary of SNCF Réseau (French Railway Network) (see [7]). These data correspond to a cruise-like flight mode in a straight line and constant ground altitude. This is the flight mode of interest for the study in this paper. For this flight mode, R and $\gamma$ should be zero, which is not the case in practice.

For simplicity reason, a bi-variate normal distribution has been fitted on these experimental data, after rejection of the outliers (see [7]). Its mean is given by $\mu ={\left[\begin{array}{cc}{\mu }_R& {\mu }_{\gamma }\end{array}\right]}^T$ with ${\mu }_R=7.47\times {10}^{-5}$ rad/s and ${\mu }_{\gamma }=1.03\times {10}^{-1}$ deg and its covariance matrix by:

$$\Sigma =\left[\begin{array}{cc}7.40\times {10}^{-4}& 1.75\times {10}^{-3}\\ 1.75\times {10}^{-3}& 8.53\times {10}^{-1}\end{array}\right]$$

(7)

This distribution will be used to sample points for $R_0$ and ${\gamma }_0$.

3.2. Uncertainties on Control Surface Deflections

As previously mentioned, a main engine failure is considered in this paper. A constant zero thrust command is therefore assumed for the descent trajectory simulation, that is ${\delta }_T\left(t\right)={\delta }_{T_0}=0$, $\forall t\in [t_0,t_f]$.

For the deflection of the control surfaces (elevators, ailerons and rudder), it is also assumed that they remain stuck in their trim position $({\delta }_{e_0},{\delta }_{a_0},{\delta }_{r_0})$ during the descent trajectory. Some noise $\Delta {\delta }_{i_0},i\in \{e,a,r\}$, is nevertheless added to these trim values, since in practice, a flapping behavior of these control surfaces has been observed on the UAV. It is defined as a zero-mean Gaussian noise of variance ${\sigma }_i^2={\rho }_i/30$, where ${\rho }_i$ stands for the amplitude range of the control surface i. The coefficient 30 has been arbitrarily chosen, but to define a variance small enough to make the new initial condition $({\chi }_{d_0},u_0+\Delta u_0)$ not to deviate too much from the computed trim point $({\chi }_{d_0},u_0)$, which is an equilibrium point for the UAV dynamics. The notation $\Delta u_0$ is used to define the noise vector ${\left[\begin{array}{cccc}\Delta {\delta }_{e_0}& \Delta {\delta }_{a_0}& \Delta {\delta }_{r_0}& 0\end{array}\right]}^{\top }$.

3.3. Monte Carlo Simulations

Let us bring in the same vector U, the 5 uncertain variables $R_0$, ${\gamma }_0$ and $\Delta {\delta }_{i_0},i\in \{e,a,r\}$, with joint density $f:{\mathbb{R}}^5\to {\mathbb{R}}^{+}$ with respect to Lebesgue measure. The computation of the impact points is then done with the deterministic process described in Section 2.1 synthesized by a scalar continuous function $\mathcal{M}:{\mathbb{R}}^5\to {\mathbb{R}}^2$. The impact position vector Z is such that $Z=[x,y]=\mathcal{M}\left(U\right)$. As U is a random vector, Z is also a random vector of unknown density $g:{\mathbb{R}}^2\to {\mathbb{R}}^{+}$ with respect to the Lebesgue measure. If we consider N independent and identically distributed (iid) samples $U_i$, $i=1..N$, with density f, we can generate N iid samples $Z_i$ of density g thanks to $\mathcal{M}$. Figure 3 shows, for instance, 2000 iid samples of $Z_i$ depending on the tuning of the wind in the function $\mathcal{M}$.

4. Density Minimum Volume Set Estimation for the Analysis of Ground Impact Footprints

A volume set is a mathematical tool that enables one to analyse the density of drone ground impacts. In this section, we describe how to define a multidimensional density minimum volume set and how to estimate them in practice, with multiple importance sampling to focus on rare events.

4.1. Definition of a Density Minimum Volume Set

The t-level set $\mathbf{L}\left(t\right)$ of the multivariate probability density g of Z is defined as follows:

$$\mathbf{L}\left(t\right)=\left\{z\in {\mathbb{R}}^2:g\left(z\right)\ge t\right\}$$

(8)

for $t\ge 0$. The level sets of the probability density function (PDF) g are defined as the mapping $ϵ$:

$$\begin{array}{ccc}\hfill ϵ:[0,supg]& \to & [0,1]\hfill \\ \hfill t& \to & {\int }_{\mathbf{L}\left(t\right)}g\left(z\right)dz=P\left(Z\in \mathbf{L}\left(t\right)\right)=P\left(g\left(Z\right)\ge t\right)=\alpha \hfill \end{array}$$

The t-level set $\mathbf{L}\left(t\right)$ of density g is the minimum volume set of probability $\alpha$ under regularity conditions [15]. A density level set can be viewed in fact as a multidimensional $\alpha$-quantile estimation.

$$V\left(t\right)=\underset{A\in {\mathbb{R}}^r}{inf}\left\{\lambda \left(A\right):P\left(A\right)\ge \alpha \right\},\alpha \in [0,1]$$

(9)

where A is a subset of ${\mathbb{R}}^r$ and $\lambda$ is a real-valued function defined on A. If $\lambda$ is the Lebesgue measure, V is a minimum volume set of probability $\alpha$.

4.2. Statistical Estimation of a Density Minimum Volume Set with MIS

In this article, we want to estimate the t-level set $\mathbf{L}\left(t\right)$ of density g for a given probability $\alpha$. The estimation principle is based on the following steps:

• Propose an estimate $\widehat{g}$ of g from a given set of samples $(Z_1=\mathcal{M}\left(U_1\right)$,...,$Z_N=\mathcal{M}\left(U_N\right)$.
• Estimate the threshold $\widehat{t}={ϵ}^{-1}\left(\alpha \right)$ with a simple binary search and determine the level set

  $$\mathbf{L}\left(\widehat{t}\right)=\left\{z\in {\mathbb{R}}^2:\widehat{g}\left(z\right)\ge \widehat{t}\right\}$$

  (10)

This estimator $\mathbf{L}\left(\widehat{t}\right)$ is a plug-in estimator of a minimum volume set [16]. To apply this 2-step procedure, it is necessary first to estimate the unknown density g. This can be done with classical Monte Carlo from samples $Z_i$, $i=1,\cdots ,N$, distributed with the unknown density g, but also with importance sampling. The principle of importance sampling is to modify the sampling distribution of $Z_i$, in order to improve the accuracy of the estimation of g on some part of its support. A comparison between Monte Carlo and classical importance sampling estimates of g is indeed performed in [17]. Depending on the value of $\alpha$, a trade-off should be made. For this purpose, we consider in this article multiple importance sampling [18] that behaves well in the heart of the distribution g, because Monte Carlo and importance sampling samples can be taken into account in the estimation of the density g. Moreover, the estimation of $\widehat{t}$ with binary search is often intractable and cannot be applied in practice, since it requires the estimation of integrals over a multidimensional domain. To avoid this difficulty, one also considers another plug-in estimator of t described in [19], based on density quantile.

To estimate the density level set $\mathbf{L}\left(\widehat{t}\right)$ with multiple importance sampling, the following computational steps are considered in this article:

• Generate a set of N independent and identically (iid) distributed samples $(U_{11},...,U_{1N})$ of density f, and apply the function $\mathcal{M}$ on these samples to determine a set of samples $(Z_{11}=\mathcal{M}\left(U_{11}\right),...,Z_{1N}=\mathcal{M}\left(U_{1N}\right))$.
• Estimate the output density ${\widehat{g}}^1$ from the samples $(Z_{11},...,Z_{1N})$ with multivariate kernel density estimate [20].
• Estimate the density h of the samples $\left\{U_{1i}\right|{\widehat{g}}^1\left(Z_{1i}\right)<\gamma \}$ for $i=1,\cdots ,N$ where $\gamma$ is set by the user.
• Generate a set of N iid samples $(U_{21},...,U_{2N})$ from density h, and applies the function $\mathcal{M}$ on these samples to determine a set of samples $(Z_{21},...,Z_{2N})$ [20].
• Estimate the density ${\widehat{g}}^{MIS}$ from the samples $(Z_{11},...,Z_{1N})$ and $(Z_{21},...,Z_{2N})$ with weighted multivariate kernel density estimate. The weight of each $Z_{ij}$ is $\frac{f\left(U_{ij}\right)}{\frac{1}{2}(f(U_{ij}+h\left(U_{ij}\right))}$.
• Estimate the threshold ${\widehat{t}}^{MIS}$ as the $(1-\alpha )$-quantile of the weighted samples $(Z_{11},...,Z_{1N}$, $Z_{21},...,Z_{2N})$.
• The level set with MIS is then estimated with

  $${\widehat{\mathbf{L}}}^{MIS}\left(t\right)=\left\{\tilde{x}\in {\mathbb{R}}^2:{\widehat{g}}^{MIS}\left(z\right)\ge {\widehat{t}}^{MIS}\right\}$$

  (11)

The choice of $\gamma$ can be made with a quantile of the samples ${\widehat{g}}^1\left(Z_{1i}\right),...,{\widehat{g}}^1\left(Z_{1N}\right)$. In this article, $\gamma$ is quantile of level $0.1$ as, from our experience, it corresponds to a good trade-off between Monte Carlo samples and extreme samples.

4.3. Application to Drone Ground Impact

An MIS algorithm for density minimum level set has been applied to the estimation of drone ground impacts with $N=1000$. In Figure 4, we present the estimation of minimum density volume set for different probabilities $\alpha$ with MIS. Importance sampling with density h has consequently increased the frequency of impacts with a high distance from the aim without requiring a large number of simulations, and thus extreme level sets are more accurate. A similar analysis is also performed in Figure 5, when a wind of 5 m·s ${}^{-1}$ with angle $-{90}^{\circ }$ is considered in the drone ground impact simulations.

5. Reliability-Oriented Sensitivity Analysis

5.1. Definition of ROSA Sensitivity Indices

Reliability-oriented sensitivity analysis (ROSA) differs from classical sensitivity analysis in the nature of the output quantity of interest under study. Indeed, sensitivity analysis focuses on the model output, whereas ROSA analyses the impact of the input uncertainty on a reliability measure. Two kinds of ROSA indices can be computed in practice [21]:

• first, target sensitivity analysis evaluates the impact of inputs over a function of the output, typically the indicator function of a critical domain. In the drone impact application, it answers the question: which uncertain inputs lead to extreme drone impact?
• second, conditional sensitivity analysis, which aims at studying the impact of inputs exclusively within the critical domain, namely, conditionally to the failure event. In the drone application, this indice determines, conditionally to an extreme impact, which uncertain inputs are the most influential.

In this article, we consider two recent target and conditional ROSA moment-independent indices ${\overline{\eta }}_i$ and ${\delta }_i^f$ [22] to analyse the influence of the $i^{th}$ component of U, $U^{\left(i\right)}$, on the scalar output quantity $\tilde{Z}={\left|\right|Z\left|\right|}_2$ for a given failure event. We propose to define here the failure event as $Z\notin \mathbf{L}\left(t\right)$, that is, the ground impact is outside a given volume set and is thus an extreme impact. The two ROSA indices are defined by the following equations for the proposed drone fallout test case with $i=1,\cdots ,5$ as there are 5 random inputs:

$${\overline{\eta }}_i=\frac{1}{2}{∥f_{U^{\left(i\right)}}-f_{U^{\left(i\right)}\mid Z\notin \mathbf{L}\left(t\right)}∥}_{L^1\left(\mathbb{R}\right)}.$$

(12)

and

$${\delta }_i^f=\frac{1}{2}{∥f_{(U^{\left(i\right)},\tilde{Z})\mid Z\notin \mathbf{L}\left(t\right)}-f_{U^{\left(i\right)}\mid Z\notin \mathbf{L}\left(t\right)}{f}_{\tilde{Z}\mid Z\notin \mathbf{L}\left(t\right)}∥}_{L^1\left({\mathbb{R}}^2\right)}=\frac{1}{2}{∥c_i-1∥}_{L^1\left({\mathbb{R}}^2\right)}.$$

(13)

where $f_{U^{\left(i\right)}}$ is the density of $U^{\left(i\right)}$, $f_{U^{\left(i\right)}\mid Z\notin \mathbf{L}\left(t\right)}$ is the density of $U^{\left(i\right)}$ conditionally to $Z\notin \mathbf{L}\left(t\right)$, $f_{\tilde{Z}\mid Z\notin \mathbf{L}\left(t\right)}$ is the density of $\tilde{Z}$ conditionally to $Z\notin \mathbf{L}\left(t\right)$, $f_{(U^{\left(i\right)},\tilde{Z})\mid Z\notin \mathbf{L}\left(t\right)}$ is the density of the couple $(U^{\left(i\right)},\tilde{Z})$ conditionally to $Z\notin \mathbf{L}\left(t\right)$ and finally $c_i$ the copula density $(U^{\left(i\right)},\tilde{Z})$ conditionally to $Z\notin \mathbf{L}\left(t\right)$. The ROSA indices ${\overline{\eta }}_i$ and ${\delta }_i^f$ take values in $[0,1]$ where the low values of these indices mean this $i^{th}$ component of U is not influential on the failure event analysis and conversely. The computation of these indices can be done with the samples generated for density level set estimation (see Section 4) and thus requires no additional calls to $\mathcal{M}$. Moreover, this methodology can be applied even if the random inputs U are dependent contrary to ROSA variance based indices [23].

5.2. Statistical Estimation of ROSA Sensitivity Indices

To practically estimate the ROSA indices ${\overline{\eta }}_i$ and ${\delta }_i^f$, the following steps are required [22]:

• Obtain $(V_1,\dots ,V_n)$ approximately i.i.d. from $f_{U\mid Z\notin \mathbf{L}\left(t\right)}$ and their corresponding value $Z_k=\mathcal{M}\left(V_k\right)$ by $\mathcal{M}$. From $Z_k$, the value of ${\tilde{Z}}_k$ is then easily computed with ${\tilde{Z}}_k=\left|\right|Z_k{\left|\right|}_2$.
• Use the sample $((V_k^{\left(i\right)},{\tilde{Z}}_k),k=1,\dots ,n,i=1,\dots 5)$ where $V_k^{\left(i\right)}$ is the $i^{th}$ component of $V_k$, to get estimates ${\widehat{f}}_{U^{\left(i\right)}\mid Z\notin \mathbf{L}\left(t\right)}$ and ${\widehat{c}}_i$ of the density $f_{U^{\left(i\right)}\mid Z\notin \mathbf{L}\left(t\right)}$ and of the copula density $c_i$ respectively. In this article, they are both estimated with the non-parametric method [20,24], but any other efficient density and copula estimation techniques can be chosen.
• Use the estimates ${\widehat{f}}_{U^{\left(i\right)}\mid Z\notin \mathbf{L}\left(t\right)}$ and ${\widehat{c}}_i$ to compute ${\overline{\eta }}_i$ and ${\delta }_i^f$ as follows:
  • for ${\overline{\eta }}_i$, estimate the one-dimensional integral ${\parallel f_{U^{\left(i\right)}\mid Z\notin \mathbf{L}\left(t\right)}-f_{U^{\left(i\right)}}\parallel }_{L^1\left(\mathbb{R}\right)}$ either by direct numerical approximation, or if $f_{U^{\left(i\right)}}$ can be sampled from, by Monte Carlo method via

    $${\widehat{\overline{\eta }}}_i=\frac{1}{N^{\prime }}\sum _{k=1}^{N^{\prime }}\left|\frac{{\widehat{f}}_{U^{\left(i\right)}\mid Z\notin \mathbf{L}\left(t\right)}\left(U_k^{\left(i\right)}\right)}{f_{U^{\left(i\right)}}\left(U_k^{\left(i\right)}\right)}-1\right|$$

    (14)

    where the $U_k^{\left(i\right)}$ are i.i.d. with common distribution $f_{U^{\left(i\right)}}$;
  • for ${\delta }_i^f$, generate $((H_{1k},H_{2k}),k=1,\dots ,N^{\prime })$ i.i.d. uniformly distributed on ${[0,1]}^2$ and estimate ${\delta }_i^f$ by

    $${\widehat{\delta }}_i^f={\displaystyle \frac{1}{2N^{\prime }}}\sum _{k=1}^{N^{\prime }}|\widehat{c_i}(H_{1k},H_{2k})-1|.$$

    (15)

The estimates ${\widehat{\overline{\eta }}}_i$${\widehat{\delta }}_i^f$ can be computed for different failure events $Z\notin \mathbf{L}\left(t\right)$ for different values of $t=t_{\alpha }$ that correspond to several minimum volume sets of probability $\alpha$. $N^{\prime }$ can be taken as large as possible, as it does not imply any calls to $\mathcal{M}$ and thus we set to $N^{\prime }={10}^4$.

5.3. Application to Drone Ground Impact Sensitivity Analysis

The algorithm proposed in the previous section has been applied with MIS samples and thus without any supplementary calls to $\mathcal{M}$ to determine the influence on the reachability of extreme drone impacts of the different components of the random vector U. The ROSA indices are computed in

Table 1. ROSA indices for ground impact analysis. Bold numbers correspond to values greater than 0.1.

(a) No Wind
ROSA indices	$t_{0.5}$	$t_{0.8}$	$t_{0.99}$
${\widehat{\overline{\eta }}}_1$	0.07	0.06	0.06
${\widehat{\delta }}_1^f$	0.04	0.04	0.03
${\widehat{\overline{\eta }}}_2$	0.05	0.03	0.03
${\widehat{\delta }}_2^f$	0.04	0.03	0.03
${\widehat{\overline{\eta }}}_3$	0.34	0.31	0.11
${\widehat{\delta }}_3^f$	0.07	0.07	0.15
${\widehat{\overline{\eta }}}_4$	0.41	0.17	0.06
${\widehat{\delta }}_4^f$	0.16	0.18	0.16
${\widehat{\overline{\eta }}}_5$	0.06	0.04	0.04
${\widehat{\delta }}_5^f$	0.06	0.04	0.04
(b) Wind of 5 m·s ${}^{-1}$ with angle $-{90}^{\circ }.$
ROSA indices	$t_{0.5}$	$t_{0.8}$	$t_{0.99}$
${\widehat{\overline{\eta }}}_1$	0.06	0.07	0.07
${\widehat{\delta }}_1^f$	0.05	0.03	0.03
${\widehat{\overline{\eta }}}_2$	0.13	0.08	0.06
${\widehat{\delta }}_2^f$	0.02	0.04	0.03
${\widehat{\overline{\eta }}}_3$	0.40	0.24	0.05
${\widehat{\delta }}_3^f$	0.16	0.22	0.23
${\widehat{\overline{\eta }}}_4$	0.34	0.13	0.06
${\widehat{\delta }}_4^f$	0.30	0.24	0.21
${\widehat{\overline{\eta }}}_5$	0.10	0.06	0.05
${\widehat{\delta }}_5^f$	0.04	0.03	0.03

for three different level sets of probability $\alpha =0.5,0.8,0.99$. The most influential variables are the third and fourth components of U, that is, the noise uncertainty on the UAV elevators and ailerons. The positions of the drone ground impact are less sensitive to the other uncertain simulation parameters in the heart of the impact position distribution. Nevertheless, when we consider more extreme impacts ($t_{0.99}$), these observations have to be mitigated. A parameter alone does not explain an extreme fallout, as the ROSA indices decrease for $U^{\left(3\right)}$ and $U^{\left(4\right)}$. A combination of parameters leads to an extreme fallout. The comparison between ${\widehat{\overline{\eta }}}_i$${\widehat{\delta }}_i^f$ here is not really relevant and does not provide much information. The sensitivity analysis gives similar results when wind is taken into account in the simulation.

6. Conclusions

The generation of extreme ground impact footprints map has been addressed in this paper for fixed-wing UAVs failure. In the proposed approach, the computation of impact points is based on simulation of a full dynamic model, including aerodynamics of the UAV and wind effect. Uncertainties accounted for in these simulations have been characterized, based on some real flight data. Monte Carlo simulations have been performed to generate footprints; however, it is not satisfying when we focus on extreme ground footprints. For this purpose, we have presented a rare-event simulation technique called multiple importance sampling to answer the issue of extreme drone ground impacts. We also show that at low computational cost, it is also possible to derive sensitivity indices to interpret the cause of extreme impacts.

Future work will include these characterizations of extreme drone impacts for the risk analysis of UAV missions. Sensitivity and robustness analysis with regard to uncertainties on some parameters of the UAV (aerodynamic coefficients) will also be considered.