# Autonomous Robotic Mapping of Fragile Geologic Features

Zhiang Chen, J Ramón Arrowsmith, Jnaneshwar Das Authors are affiliated with
the School of Earth and Space Exploration, Arizona State University, 781
Terrace Mall, Tempe, AZ 85287, USA

###### Abstract

Robotic mapping is useful in scientific applications that involve surveying
unstructured environments. This paper presents a target-oriented mapping
system for sparsely distributed geologic surface features, such as
precariously balanced rocks (PBRs), whose geometric fragility parameters can
provide valuable information on earthquake shaking history and landscape
development for a region. With this geomorphology problem as the test domain,
we demonstrate a pipeline for detecting, localizing, and precisely mapping
fragile geologic features distributed on a landscape. To do so, we first carry
out a lawn-mower search pattern in the survey region from a high elevation
using an Unpiloted Aerial Vehicle (UAV). Once a potential PBR target is
detected by a deep neural network, we track the bounding box in the image
frames using a real-time tracking algorithm. The location and occupancy of the
target in world coordinates are estimated using a sampling-based filtering
algorithm, where a set of 3D points are re-sampled after weighting by the
tracked bounding boxes from different camera perspectives. The converged 3D
points provide a prior on 3D bounding shape of a target, which is used for UAV
path planning to closely and completely map the target with Simultaneous
Localization and Mapping (SLAM). After target mapping, the UAV resumes the
lawn-mower search pattern to find the next target. We introduce techniques to
make the target mapping robust to false positive and missing detection from
the neural network. Our target-oriented mapping system has the advantages of
reducing map storage and emphasizing complete visible surface features on
specified targets.

##  I Introduction

This work is motivated by the challenges of finding and mapping fragile
geologic features such as precariously balanced rocks (PBRs) (Fig. 1 (top
left)). PBRs are fragile to being easily toppled by strong earthquake shaking,
so their existence can provide valuable information about earthquake history
(or lack thereof) in a given region [1, 2]. While current studies of such
fragile geologic features have focused on single PBR models [3], large
datasets of spatially explicit PBR fragility with detailed geomorphic and
geologic context provide a valuable assessment of sensitivity to ground
motions. Instead of having a few fragilities, it should be possible to have a
fragility histogram for the PBRs at a site. However, model assessment in a
broad statistical sense is challenging because PRBs are sparsely distributed
and difficult to inventory. Additionally, it is critical to reconstruct PBR
base since ground contact points play a significant role in characterizing
fragilities [4]. Our ultimate goal is to assess the precarious features
broadly, increasing the number of analyses and accounting for their geomorphic
setting [2, 5].

Mapping methods such as Structure from Motion [6, 7], Visual-SLAM [8, 9], and
LiDAR SLAM [10], produce dense terrain maps that include both targets of
interest and other irrelevant features. As the sparsity of targets increases,
significant amounts of computation and memory storage are spent on
reconstructing irrelevant features. Current PBR mapping work [4] bypasses this
problem by focusing on the reconstruction of single PBRs that have been found
by geologists. More importantly, mapping methods alone decouple robot
navigation from mapping, which inhibits adaptive flight paths from emphasizing
the mapping of complete visible surface features on specified targets.

Robotic search and rescue systems provide approaches to target localization
[11], however, knowing target location (i.e., geometric center) alone is
insufficient for target mapping due to ambiguity in target contours, which
might cause collision problems if the visual-motion system is not carefully
adapted. Recent object-level SLAM [12] is appealing for target mapping but
limited to room-scale environments because the target tracking and
localization depend on a depth camera, whose sensing range is inadequate for
field applications. Other object-level SLAM based on monocular cameras [13]
approximates objects with ellipsoids, which convey target pose and occupancy
information. However, to precisely map a target, additional processes are
still required. Additionally, object-level SLAM algorithms have not explicitly
considered false positive detection from neural networks. False positive
detection is less frequent when neural networks are trained from benchmark
datasets, but has commonly been seen in domain-specified, practical
applications.

Figure 1: (Top left) A precariously balanced rock (PBR) in Granite Dells,
Arizona. (Bottom left) Reconstructed world in Gazebo. (Right) Resulting mapped
PBR from the reconstructed scene. Figure 2: Perception subsystem workflow.

This paper presents an efficient, real-time, and robust system to detect,
localize, and map fragile geologic features, i.e. PBRs, in a landscape. The
autonomous system for PBR search and mapping provides access to fragility
distribution, which is valuable to earthquake studies but has not been
obtained due to the challenge of collecting large datasets. By introducing two
tracking modules, the proposed target-oriented mapping system is robust to
false positive or missing detection from a deep neural network. We implement
the system using Robot Operating System with the PX4 UAV flight
stack111https://github.com/ZhiangChen/target_mapping.

##  II Problem Description

The goal is to search for and map PBRs in a given survey region. The system
should only map the targets of interest and also emphasize complete visible
surface features of the targets including visible contact points with the
ground. The hardware platform considered within this work is a UAV equipped
with a flight controller, one/two onboard companion computers, a GPS module,
and a stereo camera. The onboard companion computers, e.g. Intel NUC i7 +
Nvidia TX2 in practice, should allow a lightweight detection neural network to
conduct real-time inference. There are bounded noises in UAV’s pose
estimation, which is feasible for a flight controller (e.g. Pixhawk with PX4
flight stack) running EKF with GPS and IMU. A stereo camera is used for target
search and mapping. The UAV’s relative elevation with respect to the terrain
during search can be greater than the depth range of the stereo camera, so the
target search only relies on one monocular camera from the stereo camera, and
the close and detailed mapping uses the stereo camera to generate dense point
clouds.

We assume there are no obstacles around PBRs during UAV mapping. The local
terrains of PBRs are within a certain slope $\Psi$ such that the local
terrains cannot present obstacles for UAV mapping. There might be false
positive or missing detections from neural networks, so the system should be
robust to such false detections. The UAV camera can capture the full body of
each PBR during search, which is feasible in practice by adjusting the
elevation of UAV search pattern.

##  III System Overview

The proposed target-oriented mapping system is composed of a perception
subsystem (Fig. 2) and a motion subsystem (Fig. 3). We start by carrying out a
lawn-mower search pattern in the survey region from a high elevation using a
UAV. Once a potential PBR target is detected by a deep neural network, we
track the target bounding box in image coordinates using a bounding box
tracking algorithm. We generate 3D points from the tracked bounding box and
resample them with our sampling-based filtering algorithm. When the 3D points
are assessed to be converging, the UAV employs a circular motion at the same
search elevation to keep filtering 3D points. A bounding cylinder (b-cylinder)
is then constructed from the converged 3D points for UAV path planning to
closely and completely map the target with SLAM. After target mapping, the UAV
resumes the lawn-mower search pattern to find the next target.

(a) Motion state diagram

(b) Controller workflow

Figure 3: Motion subsystem.

##  IV Perception Subsystem

The perception subsystem workflow is shown in Fig. 2. We deploy two tracking
modules: one for bounding boxes and one for 3D points. The bounding box
tracking module subscribes to detection messages from a neural network and
publishes tracked bounding boxes to the 3D points tracking module. The essence
of the 3D points tracking module is a sampling-based filtering algorithm where
a set of 3D points are re-sampled after weighting by the tracked bounding
boxes from different camera perspectives (see Sec. IV B). We construct a
b-cylinder from the converged 3D points and conduct close and complete target
mapping with SLAM.

###  IV-A Bounding Box Tracking

We apply the simple online and real-time tracking (SORT) algorithm [14] to
track bounding boxes from a lightweight detection neural network. When a new
detection is associated with a tracked bounding box using intersection-over-
union (IoU) and Hungarian algorithm, SORT algorithm updates the tracked
bounding box state by maintaining a Kalman filter.

When a PBR enters the image, a unique identity is registered once the PBR has
been detected by the neural network for $T_{bbox\\_hit}$ frames. This can
prevent false positive detections in a certain sense. A tracked bounding box
is deregistered if it has not been continuously detected from
$T_{bbox\\_missing}$ frames, which can alleviate missing detection issues.
Only registered bounding boxes are published to the 3D points tracker.

###  IV-B 3D Points Tracking

Having received a new tracked bounding box, the 3D points tracker first
examines if the bounding box is on the image edge. We drop the ones on the
image edge to ensure the bounding box captures the entire body of a target.
Then a set of 3D points are generated within a polyhedral cone that is back-
projected from the tracked bounding box (Fig. 4). The 3D points will be
projected to different image frames and re-sampled according to the positions
of their projected 2D points with respect to the tracked bounding boxes. The
rationale is similar to the depth filter in SVO [15], but we adapt it to
bounding boxes and use a set of 3D points to estimate a target occupancy based
on the previous work [16].

#### 1) Camera Model:

We derive the back-projection of image points, e.g. four bounding box corners,
to rays on special Euclidean group based on the models in Hartley and
Zisserman’s book [17]. A general projective camera maps a point in space
$\mathbf{L}=[x,y,z,1]^{T}$ to an image point $\mathbf{l}=[u,v,1]^{T}$
according to the mapping $\mathbf{l}=\mathbf{P}\mathbf{L}$. $\mathbf{P}$ is
the projection matrix, and $\mathbf{P}=[K|\mathbf{0}]T_{c}^{w}$, where
$K_{3\times 3}$ is the camera intrinsic matrix, and $T_{c}^{w}\in SE(3)$ is
the transformation matrix from world frame to camera frame. We can back-
project a ray from an image point,

$\displaystyle\mathbf{r}_{c}=S(\mu K^{-1}\mathbf{l})+\mathbf{b}$ (1)
$\displaystyle\mathbf{r}_{w}=T_{w}^{c}\mathbf{r}_{c}$

where $\mathbf{r}_{c}$ is a homogeneous 3D point representing a ray in camera
frame, $\mathbf{r}_{w}$ is the corresponding ray in world frame, $\mu$ is a
non-negative scalar, $S_{4\times 3}=\begin{bmatrix}I_{3\times 3}\\\
\mathbf{0}\end{bmatrix}$, $\mathbf{b}=[0,0,0,1]^{T}$, and $T_{w}^{c}\in SE(3)$
is the transformation matrix from camera frame to world frame. $S$ and
$\mathbf{b}$ convert Euclidean coordinates to Homogeneous coordinates. A ray
in world frame is represented by a set of 3D points
$\eta=\\{\mathbf{r}_{w}|\mu\geq 0\\}$. E.g. when $\mu=0$, $\mathbf{r}_{w}$ is
the position of camera origin in world frame.

Figure 4: 3D points generation and update. (Left) 3D points are generated
within a polyhedral cone that is back-projected from an enlarged bounding box
on image frame P1; 3D points are projected to a different image frame P2 and
resampled according to the positions of their 2D projection points with
respect to the new bounding box. (Right) Resampling results.

#### 2) 3D Points Generation:

We enlarge the tracked bounding boxes to compensate for errors in the UAV pose
estimation. A four-face polyhedral cone,
$\\{\eta_{w1},\eta_{w2},\eta_{w3},\eta_{w4}\\}$, is generated by four 3D rays
back-projected from four corners of the enlarged bounding boxed,
$\\{\mathbf{l}_{1},\mathbf{l}_{2},\mathbf{l}_{3},\mathbf{l}_{4}\\}$. To
randomly initialize a 3D point in the cone, we adopt the following steps.

  * •

randomly generate a direction vector from a convex combination of rays back-
projected from four enlarged bounding box corners:

$\mathbf{v}=\alpha_{1}K^{-1}\mathbf{l}_{1}+\alpha_{2}K^{-1}\mathbf{l}_{2}+\alpha_{3}K^{-1}\mathbf{l}_{3}+\alpha_{4}K^{-1}\mathbf{l}_{4}$
(2)

s.t. $\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}\sim\mathcal{U}$ and
$\alpha_{1}+\alpha_{2}+\alpha_{3}+\alpha_{4}=1$,
where $\mathcal{U}$ is a standard uniform distribution.

  * •

randomly scale the direction vector: $\mathbf{\bar{L}}=\delta\mathbf{v}$,
where $\delta$ is a random depth. $\mathbf{\bar{L}}$ is a Homogeneous 3D point
in the cone in camera frame.

  * •

convert the 3D point to world frame:
$\mathbf{L}=T_{w}^{c}(S\mathbf{\bar{L}}+\mathbf{b})$.

$\mathbf{L}$ is one randomly generated point within a cone in world frame. We
generate $m$ number of 3D points for a tracked bounding box in the matrix form
as presented in Algorithm 1.

Algorithm 1 3D Points Generation
---
|  input: $\\{\mathbf{l}_{1},\mathbf{l}_{2},\mathbf{l}_{3},\mathbf{l}_{4}\\}$,
$m$, $T_{w}^{c}$
---
output: $P_{w}$
parameter: $K$, $d_{m}$
1: $A=\mathbb{U}(4,m)$ 6: $\bar{\delta}=Diag(\delta)$
2: $B=Diag(|A|_{m})$ 7: $\tilde{P}_{c}=C\bar{A}\bar{\delta}$
3: $\bar{A}=AB^{-1}$ 8: $P_{c}=S\tilde{P}_{c}+\mathbf{b}$
4:
$C=[K^{-1}\mathbf{l}_{1},K^{-1}\mathbf{l}_{2},K^{-1}\mathbf{l}_{3},K^{-1}\mathbf{l}_{4}]$
9: $P_{w}=T_{w}^{c}P_{c}$
5: $\delta=d_{m}\mathbb{U}(1,m)$
Comments:
1). $m$ is the number of 3D points to be generated, e.g. $m=1000$.
2). $\mathbb{U}(x,y)$ is a standard uniform sampler returning a
$x\times y$ matrix.
3). $Diag(x)$ creates a diagonal matrix from a vector $x$.
4). $|\cdot|_{m}$ returns a vector of Manhattan distances of all columns.
5). $d_{m}$ is the maximum depth of 3D points along $z$ axis
in camera frame.

#### 3) Data Association:

In assigning a tracked bounding box to an existing 3D points target, each
target projects their 3D points onto the current image frame. The assignment
cost metric is the number of 2D projection points inside the bounding box. The
assignment is solved optimally using Hungarian algorithm. Additionally, a
minimum number of 2D projection points $N_{pts}$ is imposed to reject
assignment, which triggers new 3D points generation. In this study, we set
$N_{pts}=m/10$, where $m$ is the total number of 3D points for a target. If a
tracked bounding box is assigned to a target that has already been mapped, the
remaining processes of tracking and mapping are discarded.

#### 4) 3D Points Update:

Each 3D points target maintains the UAV pose at the latest time when the
target is generated or updated. Once a tracked bounding box is associated with
a target, a keyframe is selected when the difference between the current UAV
pose and the last keyframe pose maintained in the target is greater than a
threshold.

When a new keyframe is selected, we resample the 3D points of the associated
target based on the positions between the target’s 2D projection points and
the bounding box in the new image frame (Algorithm 2). The distribution
$f(\mathbf{p}|\mathbf{x})$ assesses the weights (importance) of 2D projection
points $\mathbf{p}$ on an image given a bounding box $\mathbf{x}$, and it is
formed by a combination of a Gaussian distribution and a uniform distribution:

$f(\mathbf{p}|\mathbf{x})=w_{1}h(\mathbf{p}|\mathbf{x})+w_{2}g(\mathbf{p}|\mathbf{x})$
(3)

where $h(\mathbf{p}|\mathbf{x})$ is distribution density function of a 2D
Gaussian distribution $\mathcal{N}(\mathbf{c_{x}},\Sigma_{\mathbf{x}})$, and
$g(\mathbf{p}|\mathbf{x})$ is distribution density function of a 2D uniform
distribution $\mathcal{U}$ among the bounding box $\mathbf{x}$. Parameters
$w_{1}$ and $w_{2}$ are the weights of the two distributions. The 2D Gaussian
distribution $\mathcal{N}(\mathbf{c_{x}},\Sigma_{\mathbf{x}})$ is on the
center of the bounding box, and
$\mathbf{c_{x}}=(\frac{u_{min}+u_{max}}{2},\frac{v_{min}+v_{max}}{2})$,
$\Sigma_{\mathbf{x}}=\begin{bmatrix}(\frac{u_{max}-u_{min}}{2})^{2}&0\\\
0&(\frac{v_{max}-v_{min}}{2})^{2}\end{bmatrix}$,
where $(u_{min},v_{min},u_{max},v_{max})$ are the bounding box coordinates.
The 3D points $P_{w}$ in Algorithm 2 are resampled according to their weights
from the joint distribution by importance sampling.

Algorithm 2 3D Points Update
---
|  input: $P_{w}$, $\mathbf{x}$, $T_{c}^{w}$
---
output: $P_{w}^{\prime}$
parameter: $\sigma$, $m$, $K$
1\. $\mathbf{a}=\sqrt{\sigma}\mathbb{N}(m,3)$
2\. $\bar{\mathbf{a}}=S\mathbf{a}+\mathbf{b}$
3\. $\mathbf{p}=[K|\mathbf{0}]T_{c}^{w}(P_{w}+\bar{\mathbf{a}})$
4\. $\mathbf{i}=f(\mathbf{p}|\mathbf{x})$
5\. $P_{w}^{\prime}=Importance\\_Sampling(P_{w},\mathbf{i})$
Comments:
1). $\mathbb{N}(x,y)$ is an isotropic standard multi-variate Gaussian
sampler returning a $x\times y$ matrix.
2). $\sqrt{\sigma}$ is the standard variance of update noise.
3). $f(\mathbf{p}|\mathbf{x})$ is a two-variable distribution composed of a
Gaussian distribution and a uniform distribution.

#### 5) Properties:

Two statistical metrics are used to characterize 3D points and determine a
target converging state. We use differential entropy (DE), which indicates the
compactness of a set of 3D points, to determine if a target is converging. To
compute DE, we approximate the 3D points set with a three-dimensional Gaussian
distribution, whose mean and covariance are approximated by the mean and
covariance of the 3D points. The DE is computed from the following formula for
multivariate Gaussian distribution:

$h=\frac{k}{2}+\frac{k}{2}\ln(2\pi)+\frac{1}{2}\ln(|\Sigma|)$ (4)

where $k=3$ is the dimensionality of the state vector. When $h$ is greater
than a threshold $h_{c}$, we change the target state to converging.

When a set of 3D points is updated, we want to measure the amount of
information gain from the update. One metric that captures relative entropy is
Kullback–Leibler divergence (KLD) of updated 3D points distribution with
respect to previous 3D points distribution. We compute the KLD following the
formula for multivariate Gaussian distribution:

$\displaystyle\scriptstyle
D_{KL}(\mathcal{N}_{0}||\mathcal{N}_{1})=\frac{1}{2}(tr(\Pi_{1}^{-1}\Pi_{0})+(\tau_{1}-\tau_{0})^{T}\Pi_{1}^{-1}(\tau_{1}-\tau_{0})-k+\ln(\frac{|\Pi_{1}|}{|\Pi_{0}|}))$
(5)

where $k=3$ is still the dimensionality of the state vector,
$\mathcal{N}_{0}(\tau_{0},\Pi_{0})$ is the estimated Gaussian distribution of
the updated 3D points, $\mathcal{N}_{1}(\tau_{1},\Pi_{1})$ is the estimated
Gaussian distribution of the previous 3D points, and $tr(\cdot)$ is a matrix
trace function. Because we have no prior knowledge of targets, KLD is useful
to determine if a target is converged by measuring the similarity of two 3D
points distributions. In practice, KLD of a target has to be smaller than a
threshold for $N_{KLD}$ continuous updates to be determined as converged.

#### 6) Target State Update:

The 3D points tracker processes regular tracking as described above once
receiving a tracked bounding box. When there is no message from the bounding
box tracker, the 3D points tracking process is called by a timer interrupt.
However, the timer interrupt is only in charge of target deregistration, which
is slightly different from the regular tracking process. A target is
permanently registered once mapped. Otherwise, if a target is under tracking
(the target is not converged yet), it is deregistered when it has not been
continuously updated from any tracked bounding box for $T_{pts\\_missing}$
tracking iterations including regular tracking processes or timer interrupts.
We introduce such a grace period to deal with false positive and missing
detections that are not cut out by the bounding box tracker. E.g. a 3D points
target generated from false positive detections cannot have enough updates to
get converged because the false positive detections are eliminated from some
different camera perspectives. Similarly, such grace period allows a certain
degree of missing detections as a PBR is re-detected by the neural network.

###  IV-C Target Mapping

Once a target is converged, we first examine if the target has already been
mapped by checking the distances between this target center and other mapped
target centers. If all distances are greater than a threshold, an RGBD SLAM
[9] is deployed for target mapping. When the target mapping is done, the
target state is changed to mapped, and its corresponding 3D points are
replaced by a downsampled point cloud from the mapping. The original dense
point cloud of the target is saved to a hard disk.

##  V Motion Subsystem

We start by carrying out a lawn-mower search pattern in a given survey region
from a high elevation. Once the state of a 3D points target from the 3D points
tracker is assessed as converging, the UAV deploys a circular motion at the
same lawn-mower search elevation. The purpose of this circular motion is to
keep filtering the 3D points target and ultimately estimate a b-cylinder of a
PBR from the converged 3D points. Based on the b-cylinder, the UAV conducts a
path planning to closely and completely map the target with SLAM. After target
mapping, the UAV resumes the lawn-mower search pattern to find the next
target. The motion state diagram is illustrated in Fig. 3a.

We adopt only one control system for all motion states including lawn-mower
search, b-cylinder estimation motion, and mapping motion. The control system
is illustrated in Fig. 3b. The motion states switch by altering path
generators. In the target search state, the UAV path is a hard-coded lawn-
mower pattern. For the b-cylinder estimation motion and the mapping motion,
the path generators are based on view planning methods, which are introduced
in the following subsections. The trajectory generator takes keypoints from
the path generator and generates a trajectory that minimizes the snap of 10th-
degree polynomials [18]. The geometric controller takes a desired robot’s
position and yaw heading and outputs desired robot’s body rate and thrusts
[19, 20]. The attitude controller customized with the UAV hardware
configurations by the flight controller produces desired motor speeds for UAV
propellers.

Because the geometric controller is only concerned with desired position and
yaw, we add maximum acceleration and maximum velocity constraints in the
trajectory generator such that UAV’s pose can be approximated to hovering.
That is, UAV’s roll and pitch angles are neglected in the view plannings.
Additionally, because the distance $d_{1}$ between a PBR center $C_{e}$ and
the camera is much larger than the distance $d_{2}$ between the camera and the
UAV ($d_{1}:d_{2}>10:1$), we do not differentiate camera position or UAV
position in the view plannings.

###  V-A Bounding Cylinder Estimation Motion

Because camera perspectives from lawn-mower search alone cannot always
precisely estimate PBR’s occupancy, b-cylinder estimation motion is activated
when a 3D points target is assessed as converging during lawn-mower search.
The path generator then produces a set of camera perspectives to reduce
uncertainty in the target occupancy.

We approach the view planning problem by first introducing geometric
constraints and then finding a feasible solution. Given the angle $\gamma$
between the camera $z$ axis and the horizontal plane, there are four variables
in UAV’s pose to be determined, $(x,y,z,yaw)$. We first determine $z$ by
constraining the UAV positions on the lawn-mower search plane, which is
collision-free and also parallel to $XY$ plane in world frame. Because the
view planning depends on 3D points, which might have large uncertainty before
being refined by the b-cylinder estimation motion, setting the camera $z$ axis
to pass through the geometric center of 3D points $C_{e}$ can keep the target
projection in the center of image frame as much as possible. However, this is
problematic when the camera is mounted vertically straight down. In practice,
we relax this constraint and set the camera $z$ axis to pass through a
vertical line where the 3D points center $C_{e}$ lies such that the angle
between $C_{e}$ and horizontal plane is $\gamma_{0}$, i.e.
$\gamma_{0}=45\degree$ (Fig. 5a). The UAV’s heading $yaw$ is correspondingly
determined. From the above geometric constraints, a circular path is
constructed as the path of the b-cylinder estimation motion, where the
remaining $(x,y)$ are also determined. From the current UAV pose where a 3D
points target is assessed as converging, the next keypoint from the path
generator is the closest point on the circle $\mathbf{C}_{l}$. Then the UAV
follows the circle curve $\mathbf{C}_{l}$ to keep filtering 3D points.

When 3D points are assessed as converged, we fit the 3D points with a vertical
cylinder, which is the intersection of two other vertical cylinders
$H=H_{1}\cap H_{2}$. $H_{1}$ is the smallest vertical cylinder that includes
all 3D points and has the axis traversing the center of the 3D points $C_{e}$.
$H_{2}$ has the cylinder center at $C_{e}$; its height is the six sigma value
of the 3D points along $z$ axis; its radius is the largest three sigma value
along $x$ or $y$ axis. The b-cylinder is used for view planning of target
mapping.

(a) View planning for b-cylinder estimation

(b) View planning for target mapping

Figure 5: (a) Large ellipse on the top is the circular path. $C_{e}$ is the
geometric center of the 3D points. $\gamma_{0}$ is the angle between
$\overline{C_{e}P_{t}}$ and horizontal plane. $P_{t}$ is the start position on
the circular path. (b) Three illustrated circles are the generated path for
mapping in this case. Small yellow dots represent converged 3D points of a
target. Blue cylinder is b-cylinder from the 3D points. $C_{e}$ is the center
of the 3D points (mean Euclidean distance). $r_{e}$ is the radius of the
cylinder. $r_{m}$ is the distance to be kept between the UAV and the cylinder
surface. $\beta$ is the vertical scanning field of view. $\gamma_{n}$ is the
angle between upper scanning ray and horizon plane.

###  V-B Mapping Motion

The path generator for target mapping should produce enough camera
perspectives to densely map the complete visible body of a target. We deploy a
circular path for the mapping (Fig. 5b). The circles lie on horizontal planes
that are parallel to $XY$ plane in world frame. The centers of the circles are
on the axis of the b-cylinder. We set the circle radius $r_{c}=r_{e}+r_{m}$ to
keep some distance between the UAV and the b-cylinder. The spacing and number
of the circles are determined by $r_{c}$, the upper scanning angle
$\gamma_{1}$, and vertical scanning field of view $\beta$ (Fig. 5) such that
the cylinder surface excluding the bottom will be fully covered by the camera
scanning views. We have $\beta$ smaller than the camera vertical field of view
so that there are vertical overlaps between circular scans. The first circle
is placed at the bottom where the cylinder bottom can be just covered by the
lower scanning ray. The other circles are repeatedly arranged upward until the
center of cylinder top is covered by the upper scanning ray.

###  V-C Motion State Update

As the motion state diagram shown in Fig. 3a, if 3D points have not been
converged after one circular movement of b-cylinder estimation, the system
fails to verify the 3D points target and thus deregisters the 3D points target
and resumes the lawn-mower search. When a target gets deregistered during the
b-cylinder estimation motion due to missing detections, the motion state also
switches back to the lawn-mower search.

Because the transition of the motion states is triggered by target states from
3D points tracking, when there are multiple targets with different states, we
rely on a priority table to decide the motion state. The priorities descend
from resuming lawn-mower search, to mapping motion, to b-cylinder estimation
motion. When a UAV is already on b-cylinder estimation motion state for a
target, this target has higher priority than other targets that are also
assessed as converging. In other scenarios, when there are multiple targets
with the same states, the targets are processed as a queue data structure.
Additionally, we hash target states for each target for the efficiency and
convenience of inquiry.

(a) UAV trace of Experiment I.

(b) UAV trace of Experiment II.

Figure 6: UAV traces. (a) $\overline{ABC}$ is lawn-mower search trace.
$\overline{CD}$ is b-cylinder estimation motion trace, and target is assessed
as converging at $C$ and as converged at $D$. $\overline{EF}$ is the target
mapping trace. $\overline{GH}$ is resumed lawn-mower search trace.

(a) Gazebo world for Experiment II

(b) Mapping results

Figure 7: Gazebo world and mapping results from Experiment II.

##  VI Experiment

We evaluate the target-oriented mapping system with two experiments in Gazebo
simulation using a gaming laptop, Dell G7: Intel Core i7-8750H, 16GB RAM, and
Nvidia GeForce GTX 1060 Max-Q, which has similar computing capacity with Intel
NUC i7 plus Nvidia TX2. In both experiments, we use a 3DR Iris with an RGBD
sensor ($\gamma=60\degree$) operated by PX4 SITL and MAVROS. The UAV pose
estimation is maintained by an EKF fusing data from IMU in the flight
controller and a GPS module. Sensor noises of GPS, IMU, and magnetometer are
also simulated in Gazebo. A tiny YOLO v2 is initialized from the darknet53
[21] and finely tuned with the images collected in the Gazebo worlds, and the
real-time inference is deployed using YOLO ROS [22]. We apply RTAB-MAP [9] for
RGBD mapping, and the mapping service is only activated during the mapping
motion. The vertical scanning field of view $\beta$ is $40\degree$. For the
control system, the maximum velocity is $1m/s$ and the maximum acceleration is
$1m/s^{2}$.

The first experiment is conducted to demonstrate the target-oriented mapping
pipeline for a field reconstructed from real data. We build the mesh model
with one PBR (Fig. 1 (bottom left)) using SfM with UAV imagery that was
manually collected at Granite Dells, Prescott AZ, 2019. The UAV trace with
different motion states is shown in Fig. 6a. As the mapping result presented
in Fig. 1 (right), the presented system has mapped the target and also
provided access to the complete visible surface features of the PBR including
the ground basal contact, which is crucial to calculate the PBR’s fragility
for earthquake studies.

Table I: Target-oriented Mapping System Performance

worst all detection generation converging converged mapped precision 88.7%
84.6% 88.9% 100% 100% recall 78.3% 100% 100% 87.5% 100%

The purpose of the second experiment is to examine the system with multiple
sparsely distributed targets. The 3D terrain model (Fig. 7a) is built in
Blender [23] with texture from satellite imagery and digital elevation model
from OpenTopography [24]. We import seven synthesized PBRs to the terrain
model as our mapping targets. The UAV trace is shown in Fig. 6b, and the
mapping results are displayed in Fig. 7b.

We show the system robustness by measuring the performances of key components
in Experiment II (Table. I). As false positive and missing detections are the
major cause of the problems in the following components, we inspect the worst
neural network detection performance (IoU>0.5) in one of the seven PBR
mappings, where each begins from the end of a previous PBR mapping to the
converged assessment of a current PBR. The other components are measured for
the entire seven PBR mappings based on the target states presented from the
system and the desired target states without considering mission duration. For
the false positives and false negatives of converging, we only consider the
cases that are caused by false detections from the neural network. We neglect
the cases where there are not enough detections due to the bounding box is too
close to the image edges, because they are eventually detected from other
search perspectives. This experiment shows the robustness of the system as all
seven PBRs are successfully mapped, especially with regard to the existence of
false positives and false negatives in other components.

##  VII Conclusion and Future Work

We present a target-oriented system for UAV to map fragile geologic features
such as PBRs, whose geometric fragility parameters can provide valuable
information on earthquake processes and landscape development. Such a system
provides access to autonomous PBR mapping, which was lacking in field geology
but will be essential to geologic model assessment.

Our ultimate goal is to deploy this target-oriented mapping system to actual
boulder fields and to assess the quality of 3D mapping. Additionally, 3D
semantic segmentation should be applied to extract 3D points of PBRs from
their contact terrains, which leads to PBR surface reconstruction. We have
implemented the circular motion for target mapping, which is based on an
assumption that there are no obstacles during the mapping. This can be future
improved by framing the mapping as a probabilistic exploration problem. Next-
best-view or frontier-based exploration algorithms [25, 26] can be used to
generate an obstacle-free path and accurately determine occupancy within the
bounding cylinder.

## ACKNOWLEDGEMENTS

This work was supported in part by Southern California Earthquake Center
(SCEC) award 20129, National Science Foundation award CNS-1521617, National
Aeronautics and Space Administration STTR award 19-1-T4.01-2855, and a gift
from Pacific Gas and Electric Company.

## References

  * [1] James N Brune “Precariously balanced rocks and ground-motion maps for Southern California” In _Bulletin of the Seismological Society of America_ 86.1A The Seismological Society of America, 1996, pp. 43–54
  * [2] Abdolrasool Anooshehpoor, James N Brune and Yuehua Zeng “Methodology for obtaining constraints on ground motion from precariously balanced rocks” In _Bulletin of the Seismological Society of America_ 94.1 Seismological Society of America, 2004, pp. 285–303
  * [3] David E Haddad et al. “Estimating two-dimensional static stabilities and geomorphic settings of precariously balanced rocks from unconstrained digital photographs” In _Geosphere_ 8.5 Geological Society of America, 2012, pp. 1042–1053
  * [4] Christine E Wittich, Tara C Hutchinson, J DeSanto and D Sandwell “3D reconstructions and numerical simulations of precarious rocks in Southern California”, 2018
  * [5] David E Haddad et al. “Applications of airborne and terrestrial laser scanning to paleoseismology” In _Geosphere_ 8.4 Geological Society of America, 2012, pp. 771–786
  * [6] Kendra Johnson et al. “Rapid mapping of ultrafine fault zone topography with structure from motion” In _Geosphere_ 10.5 Geological Society of America, 2014, pp. 969–986
  * [7] Zhiang Chen et al. “Geomorphological analysis using unpiloted aircraft systems, structure from motion, and deep learning” In _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2020, pp. 1276–1283
  * [8] Raul Mur-Artal and Juan D Tardós “ORB-SLAM2: An open-source slam system for monocular, stereo, and rgb-d cameras” In _IEEE Transactions on Robotics_ 33.5 IEEE, 2017, pp. 1255–1262
  * [9] Mathieu Labbé and François Michaud “RTAB-Map as an open-source lidar and visual simultaneous localization and mapping library for large-scale and long-term online operation” In _Journal of Field Robotics_ 36.2 Wiley Online Library, 2019, pp. 416–446
  * [10] David Droeschel and Sven Behnke “Efficient continuous-time SLAM for 3D lidar-based online mapping” In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 5000–5007 IEEE
  * [11] Teodor Tomic et al. “Toward a fully autonomous UAV: Research platform for indoor and outdoor urban search and rescue” In _IEEE robotics & automation magazine_ 19.3 IEEE, 2012, pp. 46–56
  * [12] John McCormac et al. “Fusion++: Volumetric object-level slam” In _2018 international conference on 3D vision (3DV)_ , 2018, pp. 32–41 IEEE
  * [13] Lachlan Nicholson, Michael Milford and Niko Sünderhauf “Quadricslam: Dual quadrics from object detections as landmarks in object-oriented slam” In _IEEE Robotics and Automation Letters_ 4.1 IEEE, 2018, pp. 1–8
  * [14] Alex Bewley et al. “Simple online and realtime tracking” In _2016 IEEE international conference on image processing (ICIP)_ , 2016, pp. 3464–3468 IEEE
  * [15] Christian Forster, Matia Pizzoli and Davide Scaramuzza “SVO: fast semi-direct monocular visual odometry” In _IEEE International Conference on Robotics and Automation (ICRA)_ , 2014
  * [16] Zhiang Chen, Sarah Bearman, J Ramon Arrowsmith and Jnaneshwar Das “Localization and mapping of sparse geologic features with unpiloted aircraft systems” In _RSS Workshop on Robots in the Wild_ , 2020
  * [17] Richard Hartley and Andrew Zisserman “Multiple view geometry in computer vision” Cambridge University Press, 2003
  * [18] Charles Richter, Adam Bry and Nicholas Roy “Polynomial trajectory planning for aggressive quadrotor flight in dense indoor environments” In _Robotics Research_ Springer, 2016, pp. 649–666
  * [19] Jaeyoung Lim “mavros_controllers - aggressive trajectory tracking using mavros for PX4 enabled vehicles”, 2019 DOI: 10.5281/zenodo.2652888
  * [20] Taeyoung Lee, Melvin Leok and N Harris McClamroch “Geometric tracking control of a quadrotor UAV on SE (3)” In _49th IEEE Conference on Decision and Control (CDC)_ , 2010, pp. 5420–5425
  * [21] Joseph Redmon and Ali Farhadi “YOLO9000: better, faster, stronger” In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 7263–7271
  * [22] Marko Bjelonic “YOLO ROS: Real-time object detection for ROS”, https://github.com/leggedrobotics/darknet_ros, 2018
  * [23] Blender Online Community “Blender - a 3D modelling and rendering package”, 2018 Blender Foundation URL: http://www.blender.org
  * [24] Giuseppe Amatulli et al. “Geomorpho90m, empirical evaluation and accuracy assessment of global high-resolution geomorphometric layers” In _Scientific Data_ 7.1 Nature Publishing Group, 2020, pp. 1–18
  * [25] Andreas Bircher et al. “Receding horizon "next-best-view" planner for 3D exploration” In _IEEE International Conference on Robotics and Automation (ICRA)_ , 2016, pp. 1462–1468
  * [26] Soohwan Song and Sungho Jo “Surface-based exploration for autonomous 3D modeling” In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 4319–4326 IEEE