Receding Horizon Navigation and Target Tracking for Aerial Detection of Transient Radioactivity

The work reported in this paper is motivated in part by application problems in the field of nuclear nonproliferation. In this context, there can be instances where one needs to quickly detect weak radiation sources that could be in transit through generally unknown, and sometimes gps-denied environments. With mavs being now a widely available and accessible technology, and with the possibility of mounting lightweight (albeit still of low efficiency) commercial off the shelf (cots) radiation detectors on them, such a detection task is now possible using aerial means.

The approach to such problems advocated here utilizes an algorithmic pipeline that combines reactive receding horizon navigation with target tracking for mavs, a visual-inertial (vi) state estimation algorithm, a ssd-MobileNetV2 based visual target identification and tracking, and a fast likelihood ratio-based binary decision making algorithm. Reactive planning is achieved through a new model predictive control (mpc)-type motion planner that fully incorporates the nonlinear vehicle dynamics of the mav (cf. Tanner and Piovesan (2010)) and utilizes real-time sensor data in the form of a point-cloud generated by an onboard rgb-d camera, to select a (probabilistically optimal) safe path toward the target within the field of view. The reactive planner solves a multi-objective optimal control problem in real time, and thus balances detection accuracy and decision-making speed against platform flight agility limitations. The product of this optimization is a dynamically compatible minimum snap trajectory that fits along the kinematic reference path. Then a nonlinear geometric controller on the mav tracks this trajectory in a receding horizon fashion.

The planning and control architecture is implemented using low-cost cots hardware, limited-range sensors, and computation devices with average capabilities. With these resources, the system has demonstrated speeds up to 4.5 m/s in cluttered indoor as well as outdoor environments. In the experimental results reported in this paper, the planner enables the mav to intercept to static or dynamic ground targets. Such interception maneuvers have been shown to increase the likelihood of correct identification of the presence of a radioactive source on the target (Sun and Tanner, 2015). In cases where the local information is insufficient to construct a feasible trajectory to the goal, the planner is designed to safely stop the vehicle. The features and capabilities of the reported reactive planning and control architecture is demonstrated in case studies of mobile radiation detection, in which mavs autonomously determine whether radiation sources are present on a ground vehicle moving in cluttered environments (Fig. 1).

A quadrotor-type mav is an inexpensive, lightweight, and agile sensor platform, suitable for many applications in areas including surveillance, aerial photography and mapping, precision agriculture, construction and defense. Although mavs with various degrees of autonomy have been deployed in these application domains, the prevailing assumption is that the environment is at least partially known, so that a motion plan can be generated a priori and then used for steering the vehicle to a desired goal. To plan its motion in unknown environments, however, a robotic vehicle needs to construct and update a (local) environment map online; recent mav literature addresses this problem using an onboard perception stack (Oleynikova et al., 2017; Hornung et al., 2013).

At a high level, the approach of this paper is different: motion planning, safe navigation, target tracking and decision-making components are integrated and co-designed to operate in concert with each other. In this regard, this paper makes technical contributions to real-time, sensor-based reactive navigation and target tracking in the form of a motion planning and control methodology applicable to problem instances involving either static or moving navigation goals. The planning and control architecture is unique because:

1. 1.

   it is applicable to cases with both static and moving goal configurations;

2. 2.

   it generates trajectories in a new way by efficiently solving a multi-objective optimal control problem;

3. 3.

   it is experimentally tested indoors and outdoors with mav speeds of up to 4.5–5 m/s; and

4. 4.

   it is demonstrated to be capable of solving a challenging radiation detection problem.

The robotics application that motivated this work was detecting weak mobile sources of radiation using aerial mobile sensors. One possibility for realizing an inexpensive fully autonomous aerial radiation detector is to utilize a quadrotor mav, mount a cots radiation counter, and allow it to maneuver itself in ways that allow it to determine as quickly as possible whether or not its given moving target is radioactive or not. While a network of static detectors can also be deployed to the detection, due to the strong (inverse square law) dependence of the process intensity on the distance between sensor and source (Nemzek et al., 2004), the static network will experience a dramatic decrease of its signal-to-noise-ratio (snr) that will significantly reduce its detection accuracy and speed.

Detecting a mobile radioactive source based on a series of discrete events associated with the arrival of particles or gamma rays (indistinguishable from each other), can be formulated as an inhomogeneous likelihood ratio test (lrt) within the Neyman-Pearson framework (Pahlajani et al., 2014). In this framework, analytical (Chernoff) bounds on error probabilities can be derived in terms of relative sensor-source distance and used as proxies for the otherwise intractable true detection error probabilities (Pahlajani et al., 2014), to facilitate the design of optimal control laws that maximize detection performance in terms of speed and accuracy (Sun and Tanner, 2015; Yadav and Tanner, 2019).

A sensor measurement event—for the case of a Gieger-Muller (gm) detector, for instance—is associated with the discrete random process where a gamma-ray (from the source or background) hits the sensitive area of the detector and records a “count.” This is mathematically modeled in a probabilistic setup (Pahlajani et al., 2014; Sun and Tanner, 2015) as follows. Assume that there is a counting process $N_{t}$ for $t\in[0,T]$, on a measurable space $(\Omega,\mathscr{F})$. In this context, $N_{t}$ essentially represents the number of times gamma-rays have hit the radiation counter located at position $x\in\mathbb{R}^{3}$, up to (and including) time $t\in[0,T]$. This counting process is assumed to be a Poisson process (Snyder, 1975). The sensor position $x$ is known over the whole interval $[0,T]$.

The detector can make two competing hypotheses, $H_{0}$ and $H_{1}$, each expressing an opinion as to whether its cumulative count rate can be attributed to the presence of a radiation source of intensity $a$ (in counts per unit time, i.e., counts per second (cps)), located at a possibly time-varying position $y\in\mathbb{R}^{3}$ which will be referred to as the target. Hypothesis $H_{0}$ asserts that the target is benign, while hypothesis $H_{1}$ claims that the target carries a source of intensity. The two hypotheses $H_{0}$ and $H_{1}$ correspond, respectively, to two distinct probability measures $\mathbb{P}_{0}$ and $\mathbb{P}_{1}$ on $(\Omega,\mathscr{F})$. With respect to measure $\mathbb{P}_{0}$, the process $N_{t}$ is a Poisson process with intensity $\beta(t)$, i.e., the intensity of naturally occurring background radiation; with respect to $\mathbb{P}_{1}$, however, the same process is Poisson with intensity $\beta(t)+\nu(t)$, where $\nu(t)$ represents the intensity of the source (whenever present) as perceived by the sensor at time $t$. The functions $\beta(t)$ and $\nu(t)$ defined on $[0,T]$ are assumed to be bounded, continuous and strictly positive (Pahlajani et al., 2014).

Function $\nu(t)$ should implicitly reflect the inverse square law dependence of the source intensity perceived by the sensor on the distance between the sensor and the source (Nemzek et al., 2004). If $\chi$ denotes the sensor’s cross-section coefficient, one possibility for expressing the detector’s perceived intensity is in the form:

A test for choosing between $H_{0}$ and $H_{1}$ faces the risk of two types of errors. One of them is a false alarm, which occurs when the outcome of the test is in favor of $H_{1}$ while $H_{0}$ is instead the correct hypothesis; the other is a missed detection in which one sides with $H_{0}$ when in fact $H_{1}$ is true. In this setting, the optimal test for deciding between $H_{0}$ and $H_{1}$ is an lrt obtained as follows (Pahlajani et al., 2014). Let $\tau_{n}$ for $n\geq 1$ denote the $n$^th jump time of $N_{t}$ (jumps occur when the sensor registers a count), and with the convention that $\prod_{n=1}^{0}(\cdot)=1$, let

be the likelihood ratio. Then for a specific fixed threshold $\gamma>0$, the test

is optimal in the (Neyman-Pearson) sense. This means that (3) minimizes the probability of a missed detection under a given upper bound constraint on the probability of false alarm. With $\mu(t)\triangleq 1+\tfrac{\nu(t)}{\beta(t)}$, constants $p\in(0,1)$ and $\eta\triangleq\log\gamma$, and the quantity

one can now extract analytical expressions for Chernoff bounds on the probabilities of false alarm $P_{F}$ and missed detection $P_{M}$ (Pahlajani et al., 2014).

If an upper limit $\alpha>0$ is set on the bound on probability of false alarm, then there exists a unique solution $p^{\ast}\in[0,1]$ for which the tightest bound on the probability of missed detection can be obtained. The exponent in the bound on the probability of false alarm and missed detection, respectively, is (Pahlajani et al., 2014)

where the dependency of $\mu$ and $\beta$ on time is suppressed for clarity, and derivative $\Lambda^{\prime}(p)=\dfrac{\partial{}\Lambda}{\partial{p}}$ is expressed as

For the optimal $p^{*}$, the $\Lambda(p^{*})$ and detection threshold $\gamma$ are related in the form $\gamma=\exp\big{(}\Lambda(p^{*})\big{)}$. Suppose now that the distance between target and sensor, $\|y-x\|$, is regulated by a control input ${u}$; then $\nu$, and consequently $\mu$, depend implicitly on ${u}$ and an optimal control problem can be formulated as follows:

Irrespective of whether $\|y-x\|$ is deterministic or stochastic, it can be shown that the optimal sensor management strategy ${u}$ for sensors is to close the distance between source and sensor as quickly as possible (Sun and Tanner, 2015; Yadav and Tanner, 2017).

Let $\mathcal{V}\in\mathbb{R}^{3}$ denote the visible space within the fov of the rgb-d sensor. This fov area contained in $\mathcal{V}$ is assumed to take the shape of a solid pyramid sector with its apex attached to the base frame of the sensor, with the depth direction of the sensor being aligned with the local $x$ (heading) frame axis of the mav. The base of the fov pyramid is outlined by the maximum range $R_{\mathrm{max}}$ of the sensor, while the side boundaries (polar and azimuth angles of the sector) are determined by the maximal viewing angles in the fov of the sensor along the local $y$ and $z$ directions. Denote $\phi_{y}$ and $\phi_{z}$ those maximum viewing angles of the sensor at its apex along the local $y$ and $z$ directions, respectively. The motion planning process also takse as input a user-specified minimum range for perception, $R_{\mathrm{min}}$. Assuming now that the apex of the fov pyramid is at local frame coordinates $(0,0,0)^{\intercal}$, any point within the fov of the sensor can be expressed in the polar coordinates as $(r,\theta,\varphi)^{\intercal}$ where $R_{\mathrm{min}}\leq r\leq R_{\mathrm{max}}$, $-\phi_{y}\leq\theta\leq\phi_{y}$ and $-\phi_{z}\leq\varphi\leq\phi_{z}$. By selecting a resolution $\delta r$ on the range and $\delta\theta$ on the viewing angles in both direction, the field of view of the sensor can be discretized and represented as an ensemble of points. Each point in this ensemble represents a potential desired next location for the mav, should a feasible path to this location exist.

The center of gravity (cog) of the mav is designated as the initial point of any candidate reference trajectory. Using the (constant) transformation between cog and sensor frame (at the fov apex), the coordinates of all the points and rays can be expressed relative to body-fixed cog frame of the vehicle.

Given the ensemble of points within the field of view, a minimum snap trajectory to each of these point is constructed. Note that the dynamics of a quadrotor mav enjoys differential flatness properties (Mellinger and Kumar, 2011), which ensure that all inputs and states can be written as functions of four (flat) outputs and their derivatives. The flat outputs are the Cartesian mav position and its yaw angle, which can be brought together in a vector $(\mathsf{x},\psi)$. Note that yaw $\psi$ is decoupled from $\mathsf{x}$ and can be steered independently.

The derivatives of $\mathsf{x}$ that are of interest are velocity $\mathsf{v}=\dot{\mathsf{x}}$, acceleration $\mathsf{a}=\ddot{\mathsf{x}}$, and jerk $\mathsf{j}=\dddot{\mathsf{x}}$. In the flat output space, the mav position dynamics can therefore be assumed to match those of a quadruple integrator with state vector $\mathsf{p}(t)=\left(\mathsf{x}^{\intercal}(t),\mathsf{v}^{\intercal}(t),\mathsf{a}^{\intercal}(t),\mathsf{j}^{\intercal}(t)\right)^{\intercal}$. The objective of the receding horizon planner now is to generate smooth trajectories $\mathsf{p}_{[N]}$ and $\mathsf{p}_{[N+1]}$, defined over the course of two consecutive planning cycles each of duration $\delta t$ and indexed $N$ and $N+1$,

which always remain within $\mathcal{F}$ and satisfy the boundary condition $\mathsf{p}_{[N]}(t+\delta t)=\mathsf{p}_{[N+1]}(t+\delta t)$, while being dynamically feasible, i.e., $f_{\mathrm{min}}\leq f\leq f_{\mathrm{max}}$ and $\|\Omega\|\leq\Omega_{\mathrm{max}}$.

Let $\mathsf{T}$ be the free trajectory terminal time (which will henseforth be referred to as the planning horizon) and denote $\mathsf{p_{0}}=\left(\mathsf{x}_{0}^{\intercal},\mathsf{v_{0}}^{\intercal},\mathsf{a_{0}}^{\intercal},\mathsf{j_{0}}^{\intercal}\right)^{\intercal}$, $\mathsf{p_{\mathsf{T}}}=\left(\mathsf{x_{T}}^{\intercal},\mathsf{v_{T}}^{\intercal},\mathsf{a_{T}}^{\intercal},\mathsf{j_{T}}^{\intercal}\right)^{\intercal}$ the trajectory boundary conditions at $t=0$ and $t=\mathsf{T}$, respectively. Then let $\mathsf{u}(t)=\dfrac{\operatorname{d\!}{{}^{4}}\mathsf{x}(t)}{\operatorname{d\!}{t^{4}}}$ denote snap, and treat it as the control input for the flat space dynamics. For a relative weight parameter $k\in\mathbb{R}_{+}$, the free terminal time optimal control problem the solution of which are the desired reference trajectories, is formulated as:

The cost function of the optimal control problem (11) is a linear combination of two performance factors: (i) the incremental cost associated with the terminal time (time optimality), captured by the constant integrand term; and (ii) the incremental cost that penalizes the cummulative magnitude of snap $\mathsf{u}$ along the trajectory. By tuning $k$, one can adjust how aggressive the reference trajectory obtained is. Smaller values for $k$ place a smaller penalty on tracking time and therefore result to slower motion. gazebo simulation data at various $k$ and speeds had been collected to fit a relation between the maximum speed along a candidate trajectory and parameter $k$. This relation comes out to be an interpolated curve of the form $\mathsf{v}_{\mathrm{candidate}}=\alpha\,k^{1/\beta}$ and has been found to work well in practical scenarios to guide the selection of the cost weight $k$ based on the maximum robot speed afforded for the mission at hand. The cost weight $k$ is particular to the candidate trajectory and varies for the different candidate trajectories in the ensemble, since the maximum speed afforded along a candidate trajectory $\mathsf{v}_{\mathrm{candidate}}$ itself varies in proportion to the ratio of the candidate trajectory’s length to the length of the longest candidate trajectory. This leads to smaller trajectories having lower top speeds, making the mav automatically slow down in the vicinity of dense obstacles.

Denoting $t$ the time elapsed since the start of the whole planned maneuver, $d$ the vehicle’s remaining distance to its goal, $r$ being the distance of the point in the ensemble from the camera focal point, $v_{\mathrm{max}}$ the desired maximum mav speed, and $k_{t}$ and $k_{d}$ being positive tuning parameters, the velocity used to calculate the weighing factor $k$ in (11) to generate a particular candidate trajectory is given by\useshortskip

Compared to alternative trapezoid velocity profiles (Mellinger and Kumar, 2011), (12) produces a velocity profile for the entire remaining quadrotor trajectory that is also effective in scenarios involving moving target interception, in which the vehicle needs to adjust its speed to match that of its target while at the vicinity of the latter. A more detailed discussion of the effect of this velocity profile on tracking dynamic targets follows in the section on target tracking.

To solve (11) one can utilize Pontryagin’s Minimum principle (Athans and Falb, 1966, Chapter 6) as follows. Let $\mathsf{\lambda_{x}}$, $\mathsf{\lambda_{v}}$, $\mathsf{\lambda_{a}}$ and $\mathsf{\lambda_{j}}$ be the costate vectors. Each such vector has three components, one for each spatial direction $x$, $y$, and $z$. Let $\langle\cdot,\cdot\rangle$ denote the standard inner product between of two vectors, and express the Hamiltonian $\mathsf{H}$ for this problem as

The optimal solution is derived from the condition $\mathsf{H}(\mathsf{x^{*}},\mathsf{u^{*}},\mathsf{\lambda^{*}},t)\leq\mathsf{H}(\mathsf{x^{*}},\mathsf{u},\mathsf{\lambda^{*}},t)$, which, since the Hamitonian is quadratic in the control inputs, leads to

The costate dynamics now is

Fortunately, this problem affords independent solution along each direction. The following equation set, where $c=(c_{x0},\cdots,c_{x7})$ denotes an $8\times 1$ vector of constant coefficients, illustrates the solution along in $x$ direction; the other two directions feature identical polynomials:

The optimal trajectory, therefore, is a 7^th order polynomial in time. Enforcing the boundary conditions at $t=0$ gives $c_{x0}=\mathsf{p}_{x0}$, $c_{x1}=\mathsf{v}_{x0}$, $c_{x2}=\mathsf{a}_{x0}$ and $c_{x3}=\mathsf{j}_{x0}$, while the remaining coefficients, $c_{x4}$ through $c_{x7}$, are obtained from the boundary conditions at $t=\mathsf{T}$:

expressing the optimal trajectory coefficients as a function of the (yet unknown) free terminal time $\mathsf{T}$.

The free terminal time $\mathsf{T}$ can be determined as follows. With the $c_{x}$ trajectory coefficients are explicitly expressed in terms of $\mathsf{T}$, one substitutes and writes the control input at the terminal time as

(control inputs at $\mathsf{T}$ in $y$ and $z$ are similarly obtained). Velocity, acceleration and jerk at time $\mathsf{T}$ are all set to zero, while the transversality condition (Athans and Falb, 1966) at $\mathsf{T}$ requires $\mathsf{H(T)}=0$. Combining these conditions, with (14) results to

This is essentially an 8^th degree polynomial equation which can now produce $\mathsf{T}$. Indeed, if for the sake of brevity we set

(work similarly for $y$ and $z$ using indices $2$ and $3$, respectively), and then substitute back we obtain the polynomial equation

This equation can be efficiently solved numerically, returning three pairs of complex conjugate roots and a single pair of real roots, of which one is negative and the other is positive—the only acceptable root.

The vehicle’s yaw angle $\psi$ is dynamically decoupled from its position. A reference trajectory for yaw can be constructed in a similar way, given that the relative degree of flat output $\psi$ is two. In the flat output space, therefore, the yaw dynamics can be expressed as as a double integrator. Let us define the yaw state vector $\Psi\triangleq(\psi,\dot{\psi})^{\intercal}$, treating $w\triangleq\ddot{\psi}$ as the virtual control input for this dynamical subsystem. For the yaw trajectory we have boundary conditions on both sides, denoted $\Psi(0)=\Psi_{0}\triangleq(\psi_{0},\dot{\psi}_{0})$ and $\Psi(\mathsf{T})=\Psi_{\mathsf{T}}\triangleq(\psi_{\mathsf{T}},\dot{\psi}_{\mathsf{T}})$. The reference yaw trajectory is obtained as a solution of

which is a 3^rd order polynomial of the form

with coefficients given as

Figure 3 shows a sample evolution of the cost functional in (11) for some particular set of boundary conditions. The cost associated with the terminal time (dotted curve) increases while the integral of the snap (dashed curve) reduces with time. The vertical line represents the positive real root $\mathsf{T}$ of (16) which minimizes the total cost (solid curve).

In the context of the particular case study considered here, the best radiation detection performance has been shown to be achieved when the robot closes the distance between its sensor and the source as fast as possible (Sun and Tanner, 2015; Yadav and Tanner, 2019). This implies that aggressive reference trajectories are preferable, motivating us to select $k$ based on the maximum speed limitations of the mav. With this in mind, the methodology outlined above can thus provide motion plans that would be (a) feasible given the dynamic constraints of the vehicle, (b) as aggressive as possible, and (c) high-performance in terms of radiation detection accuracy.

While the reference trajectory is dynamically feasible by design, its conformity to actuation constraints is verified after its generation (Mueller et al., 2015) (see Fig. 4). Here, we have actuation (upper) bounds on the magnitude of the input thrust $f=m\mathinner{\!\left\lVert\mathsf{a}-\mathsf{g}\right\rVert}$ and on the square of the sum of roll and pitch angular velocities in terms of jerk and thrust, $\nicefrac{{\mathinner{\!\left\lVert\mathsf{j}\right\rVert}^{2}}}{{f^{2}}}$. These actuation constraints stem from (a) equipment limitations, and (b) the maximum deceleration that the vehicle can undergo during an emergency stopping maneuver. In the reactive navigation architecture described here, emergency stopping maneuvers are engaged when the mav cannot find a safe path in its free workspace. Finally, a linear velocity constraint is imposed in order to reduce motion blur which would otherwise affect the robot’s visual inertial odometry (vio).

Once the ensemble of candidate reference trajectories is generated, the candidate trajectories that intersect with any obstacle is excluded from the ensemble. The methodology is motivated from earlier work of Yadav and Tanner (2020), and has been modified to suit the new trajectory generation methodology.

The obstacle point cloud $\mathcal{P}$ is first converted into a kd-Tree and each trajectory is discretized and represented in the form of a finite (size $n$) sequence of points. Thereafter a query is run to find points of the kd-tree that lies within a ball of radius $r$ from each of these points on the candidate trajectory. Parameter $r$ is chosen so that it can fully enclose the quadrotor, with a suitable—based on how conservative with respect to the practical risk collision due to uncertainty or disturbances—safety margin. A candidate trajectory is said to be intersecting with an obstacle if any point (among the n points in which it is discretized) on it has any point of kd-tree within this ball of radius $r$; such a trajectory is excluded from the ensemble.

Among the collision-free candidate trajectories, the optimal one should strike an acceptable balance between safety (against collision), and speed of convergence to the goal point. An illustration of the subsequent process of selecting this optimal reference trajectory is given in Fig. 5 (for details on the associated computational requirements, see Section 7). In lieu of a global planner which would otherwise sketch a complete path from initial to final configurations, the reactive local planner defines an intermediate point as the point in the (collision free) ensemble closest to the final goal (denoted IP in Fig. 5(a). It then assigns a cost to each trajectory, in the form of a linear weighted sum of two cost components: the first component is based on the distance of each trajectory end-point to the intermediate point, normalized over the maximum distance; the second component is a normalized collision cost that captures how close the trajectory comes to $\mathcal{P}$.

Denote the total number of collision-free trajectories and the Euclidean distance between the end point of the $i$^th trajectory and the intermediate point, $p$ and $d_{i}$, respectively. Set $d_{\mathrm{max}}\triangleq\max_{i}d_{i}$, and let $\hat{r}\geq r$ be an additional safety margin (over the radius around detected obstacles). The minimum distance $\rho_{i}$ of trajectory $i$ to obstacles is found by quering a kd-Tree over $\mathcal{P}$ and minimizing over the query response. With these data at hand, the collision cost for trajectory $i$ is found as\useshortskip

The collision cost function normalizes the cost of each trajectory into the $[0,1]$ interval. Thus any trajectory that touches a ball around its nearest obstacle is assigned a collision cost of 1, while any trajectory that lies at least $\hat{r}$-away from every obstacle incurs zero cost. All other trajectories are assigned costs within the $(0,1)$ interval. The end point of the trajectory with the lowest total cost becomes the local goal (i.e., within the fov) for the planner (marked with a green dot in Fig. 5(c)).

Selecting positive weights $k_{1},\,k_{2}\in(0,1)$, the cost contributions of trajectory $i\in\{0,\ldots,p\}$ due to its proximity to the intermediate point and obstacles are combined into an aggregate cost expression

The trajectory with the least such cost is selected as the reference trajectory between the robot’s current position and the local goal.

The mav tracks a small portion of the reference trajectory, hereafter referred to as the control horizon. The length of the control horizon is dependent on the vehicle’s speed, its computational capabilities, and its sensor update rate. Each time the mav receives a sensor update, it generates a new reference trajectory (segment) and appends it to the end of the segment it is currently tracking. By design, the transition between subsequent reference trajectory segments is smooth. The process in between sensor updates constitutes a replanning cycle. This process is illustrated in Fig. 7 which shows different trajectories generated from the starting point until the end point. The complete implementation for planning, control and state estimation is open-source.¹¹1 https://github.com/indsy123/quad_navigation_and_target_tracking

The motivating application for this work is the detection of weak radiological material on moving ground vehicles, and for this to be achieved via aerial means of sensing and measurement in a timely manner, the mavs need to get close to their target as quickly as possible. To this end, the methodology of Section 5 is extended to dynamic target interception and tracking scenarios. Necessary extensions include the capability of the mav to autonomously detect its target within its fov, and estimate its target’s relative position and speed based on visual means. For this task an open source implementation of MobileNetV2 (Sandler et al., 2018) deep neural network in TensorFlow (Abadi et al., 2016) is utilized and combined with ssd (Liu et al., 2016) for object detection.

This implementation of ssd-MobileNetV2 has been trained on the standard coco dataset. The desired neural network is then trained via transfer learning on a dataset of approximately 500 images of the Jackal mobile robot from Clearpath Robotics (Fig. 6(a)). These images were collected from an Intel RealSense D435 camera onboard the mav in different indoor and outdoor environments, under different lightening conditions and varying background. The network utilizes a $300\times 300$ fixed image resizer to increase the inference speed. Adam optimizer was used with initial learning rate of $2\times 10^{-4}$ and subsequently $10^{-4}$, $8\times 10^{-5}$ and $4\times 10^{-5}$ after 4500, 7000 and 10000 steps, respectively, for a total $2\times 10^{4}$ steps. All other parameters are kept to their default values. The network is trained on nvidia gtx 1060 gpu and the trained model is subsequently converted into a TensorRT model to run fast inference on nvidia Jetson-Nano gpu onboard the mav.

Given a bounding box obtained from the dnn, the onboard rgb-d sensor data allow direct measurement of relative position of the target with respect to the quadrotor. First, all sift features are extracted for two consecutive detection boxes and matched using a flann-based matcher utilizing a Lowe’s ratio test (Lowe, 2004), and a ransac (ransac) based outlier rejection was implemented. Utilizing the disparity image from the rgb-d sensor, the $(u,v)$ position of all the inliers (i.e. the features on the target) in the focal plane can be converted to 3d position of the feature with respect to the camera by utilizing the disparity map. The average of these 3d positions provides an aggregate measurement of relative position of the target with respect to the mav. This aggregate relative position measurement is used by a Kalman filter, which based on a constant-acceleration motion model for the target, returns 3d position estimates of the target with respect to the mav’s cog and feeds it to the motion planning algorithm. The entire target position estimation workflow is showcased in Fig. 6(b).

The effectiveness of the velocity profile prescribed in (12) is pronounced in the case of intercepting a moving target. Note that in case of a static navigation end goal, the first factor (the time-dependent erf) increases slowly from $0$ and converges to $1$, while the second term (distance-dependent erf) starts from $1$ and converges to $0$; this steers the mav so that it starts at lower speed, progresses to near maximum allowed speeds during most of the navigation maneuver, and then smoothly slow down near its goal, thus mimicking a trapezoidal velocity profile. In contrast, in the case of moving target interception the distance-dependent erf in (12) converges instead to a fixed strictly positive value, allowing the mav align its velocity with that of the target and maintain a fixed distance to it.

The reported planning and control framework was first tested in simulations based on the RotorS package (Furrer et al., 2016), within a Poisson forest-like environment with obstacle densities of 18 and 36 obstacles in a 100 m² area. A typical rviz visualization of the resulting algorithm execution is shown in Fig. 7. As the mav flies from left to right, the figure marks in green color the different trajectories generated within each planning cycle, while the fov grid point ensemble is shown in white toward the end of the vehicle’s path. Figure 8 shows the probability of mission success as a function of mav speed at two different obstacle densities (cf. Karaman and Frazzoli (2012)).

The custom-build mavs used for experimental testing shown in Fig. 9 are based on a dji Flamewheel F450 frame. Computationally, one of them features an on board Intel nuc Core i7-8650U quad core cpu@1.9 GHz$\times$4 while the other has an Intel nuc Core i5-7300U dual core cpu@2.6 GHz$\times$2. Both uses 16GB ram, 128GB ssd and a Pixhawk flight controller that is given the desired thrust magnitude and rate commands, which are then tracked using an onboard body rate controller. A point cloud is generated by an Intel RealSense-D435 depth rgb-d camera (640$\times$480 pixel, 30 Hz) while the RealSense-T265 vi-sensor (2 848$\times$800 pixel 30 Hz cameras with a 200 Hz imu) is used for inertial odometry. This lightweight sensor package provides reliable depth information for up to 5 m. A voxel filter reduces the density of the generated point cloud to a uniform size of 0.125 m, which is considered adequate for typical obstacle avoidance purposes. Open-vins provides state estimation and real time odometry at 30 Hz. Ultimately, the pipeline’s computational bottleneck is the inference script that can only be run at 15 Hz on the cpu of the intel nuc, and therefore the entire receding horizon planning and target tracking software is constrained at 15 Hz.

Over five different runs, each of overall trajectory length of 25 m, in both indoor and outdoor environments the 95% quartile of execution time is well below 0.02 seconds (Fig. 10). These execution times are almost half of those obtained on those systems in earlier studies (Yadav and Tanner, 2020) that uses an optimization based trajectory generation methodology. In this configuration, the mav flies safely among moderately dense outdoor obstacle environments at speeds of 4.5–5 m/s. These speeds surpass or are comparable to those reported in recent literature (Gao and Shen, 2016; Sikang Liu et al., 2016; Fragoso et al., 2018) without using any high-end and expensive sensors and their achievement is attributed to the ability to replan faster.

The video attachment that accompanies this paper features a number of experiments, the first two of which show the mav to navigate at speeds of 2 m/s and 4.5 m/s in an outdoor environment for a total of 40 m and 60 m respectively. The video showcases first and third person views of the cluttered scene traversed by the mav, alongside an rviz representation of generated trajectories amongst the sensed obstacles. The octomap featured in the video is provided for illustration purposes only to mark obstacle locations along the path of the vehicle, and is not used for planning. The third experiment included in the video attachment demonstrates the mav’s capability to avoid moving obstacles in its field of view, while the fourth showcases receding horizon navigation with combined obstacle and target tracking abilities. The latter utilizes a neural network with feature matching between two subsequent images.

Ultimately, the mav’s speed will be limited primarily by the computational capabilities and the need to implement a safe emergency stop in the absence of an effective motion plan—a possibility which cannot be eliminated in purely reactive and local planning methods. Since purely reactive approaches are supposed to rely exclusively on local information, convergence to the navigation goal cannot be guaranteed for all possible scenarios. (The reported reactive planning architecture can however be integrated with an exact global planner, and overcome these limitations at the price of using global information (Yadav and Tanner, 2021b).) This investigation naturally exposes the limits of purely reactive motion planning approaches. It is expected that knowledge of those limits, can guide the development of hybrid (local & global (Yadav and Tanner, 2019)) mav motion planning methodologies destined for deployment in environments where uncertainty is reasonably well characterized, in order to complement each other and operate robustly in real-world scenarios.

While on one hand the of perception from trajectory generation can induce latencies, their co-design using powerfull on board computation and the incorporation of pre-trained neural networks for trajectory generation can further boost vehicle speeds to impressive levels (Loquercio et al., 2021). To reach speeds up to 10 m/s, however, one would also require a very favorable power-to-weight ratio (e.g. $\sim$4.0 (Loquercio et al., 2021)), which may be challenging to maintain depending on the mission-mandated sensor payload. In this paper, the mavs featured a power-to-weight ratio in the 1.55–1.75 range. Besides an effort to achieve a more favorable power-to-weight ratio, we postulate that additional speed improvements can be achieved with the use of event cameras (Falanga et al., 2020).

The radiation detection algorithm is based on a Neyman-Pearson based, fixed time interval binary hypothesis test (refer to Pahlajani et al. (2014) for a more detailed exposition). At the heart of the test is a likelihood ratio of a statistic $L_{T}$ (2) calculated based on the history of relative distance between the airborne radiation sensor and the hypothesized source in addition to the aggregated counts over the sensor’s (predetermined) integration time interval $T$. This likelihood ratio is compared against a fixed threshold value $\gamma$ that also depends on the relative distance and the acceptable bound on the probability of false alarm $P_{\mathsf{FA}}$. The optimal value of $p^{*}$ is obtained by solving (5) and then the threshold is calculated by evaluating (7) at $p^{*}$. The remaining parameters for this test, including background naturally occurring (background) radiation rate and sensor characteristics, are determined by computational and experimental calibration processes.

Our detection tests involve a sequence of controlled experiments in which the efficiency of the radiation sensors and the variation of its expected count rate as a function of its distance to the source is estimated. The mavs featured in Fig. 9 were deployed in both indoor and outdoor experiments, where their task was to locate and intercept a ground target (the remotely controlled ClearPath Robotics Jackal in Fig. 1) moving along an unspecified path with unknown but bounded speed. The ground robot carried an approximately 8 $\mu$Ci radioactivity source which the mavs had to detect.

Each experiment involved the mav tracking the ground vehicle for certain time $T$. To minimize radiation exposure risk, Monte Carlo simulations using gazebo have been performed to compliment the experimental validation. The counts were generated using thinning algorithm in the simulations (Pasupathy, 2009). These simulations were used for the estimation of the minimum length $T$ of the time interval for radiation measurement collection that is more likely to yield a confident and accurate classification of the target relative to its radioactivity. This overall process suggested a sensor integration window set at $T=100$ seconds for the gm-10 counter (at the median of the distribution with 5% and 95% percentiles estimated at 71 and 136 seconds, respectively), for a radiation source of activity around $8.2\,\mu$Ci. During that fixed time interval, the distance between sensor and source should be at most within 3–2.5 m, otherwise the source emissions completely blend into the background, rendering detection extremely unlikely given the mav’s operating time constraints imposed by on-board power. For that same source and range to sensor, the gm-90 counter appeared to need approximately $T=70$ seconds, with a 5% and 95% percentiles at 65 and 96 seconds, respectively.

The receding horizon planning and control strategy of Section 5 ensures that the mav closes this distance as fast as possible, thus enabling the onboard Geiger counters to collect informative measurements. The mav has to maintain a certain minimum distance from the target to keep it in its camera’s limited fov; as a result, the relative distance should not be arbitrarily reduced. Variations in the relative distance between sensor and source can be attributed to motion perturbations, as the (remotely operated) target performs avoidance maneuvers. Although generally robust, the neural network does not furnish guarantees against false target identification, and this can contribute to relative distance estimate outliers (see e.g.  Fig. 11, around the 85^th second).

Figure 11 presents the results of one radiation detection experiment conducted in an indoor environment (an abandoned warehouse) using the mav that has GM-10 counter. It shows the evolution of the estimate of the relative distance $d$ between the mav and the ground robot as the latter moves with unknown and time-varying speed. The relative distance is estimated in real time via the target tracking pipeline described in the section on target tracking. The dashed curves in Fig. 11 indicate the evolution of Chernoff bounds on the probability of false alarm, $P_{\mathsf{FA}}$, and probability of missed detection $P_{\mathsf{M}}$ (Pahlajani et al., 2014). The bound on the probability of false alarm appears to drop below the acceptable upper limit after approximately 50 seconds from the start of the experimental run, after which the bound on the probability of missed detection $P_{\mathsf{M}}$ also starts to slowly decrease monotonically—the latter is a decreasing function of the sensor integration time and distance between sensor and source (Pahlajani et al., 2014). The graph of the logarithm of the likelihood ratio $L_{T}$ over the detection threshold $\gamma$ over time is marked in red; this process is stochastic because it depends directly on the arrival time of gamma rays on the sensor. The initial segment of the red curve corresponds to the initial time period during which the constraint on $P_{FA}$ has not been satisfied and $\log\nicefrac{{L_{T}}}{{\gamma}}$ has been kept at 0. The experiment is concluded at $95.54$ seconds and the likelihood ratio $L_{T}$ exceeds its threshold value at $89.8$ seconds indicating the presence of the radiation source on the ground target (marked with a black circle in the plot). The likelihood ratio had actually crossed the threshold before that time, but the experiment was continued because that event was observed significantly earlier than the recommended sensor integration window.

Figures 12(a) and 12(b) showcase two different runs where the mav featuring the gm-90 counter was utilized. The experimental run of Fig. 12(a) shows an instance where the mav did not have enough time to detect the source. This experiment was performed in the same indoor facility as that used for the run of Fig. 11. Here, the radiation sensor integration window is 56 seconds. The bound on the probability of missed detection is then still around 0.6, comparable to the conditions under which the detection of Fig. 11 was achieved, but this $T$ is below the 5% percentile for the recommended exposure time.

Figure 12(b) depicts the results of a longer chase by the mav carrying the gm-90 counter conducted outdoors. This time, the integration window was extended to more than 200 seconds. In addition to the effect of sensor integration window length on detection accuracy, Fig. 12(b) shows more clearly the evolution of the bounds on the decision test’s error probabilities $P_{\mathsf{FA}}$ and $P_{\mathsf{M}}$. At the time of decision, the bound on the probability of miss, $P_{\mathsf{M}}$ is almost zero, indicating very high probability for accurate decision-making. Although the statistic $\log\nicefrac{{L_{T}}}{{\gamma}}$ crosses becomes positive for the first time shortly after 70 seconds, at that time the bound $P_{\mathsf{M}}$ is around 0.3. It is of interest that towards the end of the integration window, the statistic $\log\nicefrac{{L_{T}}}{{\gamma}}$ decreases, most likely due to the target being able to open up its distance with respect to the pursuing mav—which by that time was experiencing a drop in its power reserves; same trend can be noticed in Fig. 12(a).

The accompanying video attachment includes two experiments of target tracking with radiation detection (experiments #5 and #6). These are cases where the Jackal robot is steered manually while carrying an $8.2\,\mu$Ci source. The plot at the bottom left of the video represents graphically the evolution of the detection parameters as depicted in Figs. 11 and 12(b). To reduce the size of the video file, only selected portions of these experimental runs are included and the video is accelerated at four times the normal speed.