\section{Introduction \label{sec:Introduction}}

The homogeneous electron gas (HEG) has a special place in the history
of the study of many-electron systems in general, and of
density-functional theory in particular.\citep{Hohenberg1964,
  Kohn1965} In fact, the development of accurate exchange-correlation
functionals typically begins with the local (spin) density
approximation (LDA), whose construction is based on the
exchange-correlation energy of the HEG. This is then modified by an
enhancement factor that depends on the gradient of the density in the
generalized gradient approximation (GGA); the mega-GGA approximation
adds further dependence on the local kinetic energy density and/or the
Laplacian of the electron density.\citep{Becke2014, Jones2015,
  Mardirossian2017a}

LDAs, GGAs, and meta-GGAs form the first three rungs of the so-called
Jacob's ladder of density functional theory,\cite{Perdew2001} each
rung generally leading to approximations of better accuracy.  Although
GGAs and meta-GGAs add more physical information into the density
functional approximation (DFA), they are typically constructed to
maintain exactness for the exchange-correlation energy of the HEG. In
fact, it can be even argued that this is one of the most important
exact conditions that a functional should fulfill.

In this work, we investigate the accuracy of an ansatz, which alike
the LDA is derived from considerations of the HEG only, but which adds
a further dependence on the local kinetic energy density, similarly to
meta-GGAs. These functionals, which we term meta-LDA functionals, thus
constitute a new rung on Jacob's ladder of functionals, that is shown
to fall in accuracy between LDAs and GGAs.

The work is organized as follows. We will describe the theory behind
the meta-LDA approach in \secref{Theory}, and the implementation of
the meta-LDA functionals and the details of our computations in
\secref{Computational-Details}. The accuracy of the novel functionals
is then assessed by benchmarking exchange energies of atoms and
atomization energies of molecules in \secref{Results}. A brief
summary and conclusions are presented in \secref{Conclusions}. Atomic
units are used throughout the manuscript, unless specified otherwise.

\section{Theory \label{sec:Theory}}

The LDA for the exchange energy is derived for the HEG as\citep{Bloch1929, Dirac1930}
\begin{equation}
E_{x}^{\text{LDA}}[n]=-C_{x}\int n^{4/3}(\boldsymbol{r}){\rm d}^{3}r\label{eq:lda}
\end{equation}
where
\begin{equation}
C_{x}=\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}.\label{eq:Cx}
\end{equation}
The kinetic energy density of the gas is also known,
\begin{equation}
\tau^{\text{HEG}} = C_F n^{5/3}\label{eq:tau-lda}
\end{equation}
where
\begin{equation}
C_{F}=\frac{3}{10}(3\pi^{2})^{2/3}.\label{eq:Cf}
\end{equation}
Since \eqref{tau-lda} establishes a link between the kinetic energy
density and the electron density, \citet{Ernzerhof1999a} proposed an
exchange functional similar to \eqref{lda} where \eqref{tau-lda} is
used to replace the local density dependence by
\begin{equation}
  \tilde{n}(\boldsymbol{r})=\left[\frac{\tau(\boldsymbol{r})}{C_F}\right]^{3/5}\label{eq:ntau}
\end{equation}
yielding the local $\tau$ approximation (LTA) exchange functional
\begin{align}
E_{x}^{\text{LTA}}[\tau] & =-C_{x}\int\left[\frac{\tau(\boldsymbol{r})}{C_{F}}\right]^{4/5}{\rm d}^{3}r.\label{eq:lta}
\end{align}

Based on the work of \citeauthor{Ernzerhof1999a}, \citet{Eich2014}
suggested another exchange functional where only the energy per unit
particle is written as a function of the fictitious density of
\eqref{ntau}, yielding the tLDA exchange functional
\begin{align}
E_{x}^{\text{tLDA}}[n,\tau] & =-C_{x}\int n(\boldsymbol{r})\tilde{n}^{1/3}(\boldsymbol{r}){\rm d}^{3}r.\label{eq:tlda}
\end{align}

In this work, we show the power of this idea by generalizing the
approach of Ernzerhof, Scuseria, Eich, and Hellgren. We thus replace
the electron density by an effective density $n^\text{eff}({\bf r})$
formed as a weighted combination of the electron density
$n(\boldsymbol{r})$ and the fictitious density computed from
$\tau(\boldsymbol{r})$ as
\begin{equation}
n(\boldsymbol{r})\to
n^{\text{eff}}(\boldsymbol{r})=\tilde{n}^{x}(\boldsymbol{r})n^{1-x}(\boldsymbol{r}).\label{eq:wdens}
\end{equation}
This form interpolates between the LDA ($x=0$), tLDA ($x=1/4$) and LTA
($x=1$) in the case of the exchange functional. Furthermore, it can
also be employed within any LDA correlation functional, allowing us to
generate a complete exchange-correlation ansatz.

We note here that the family of functionals generated by \eqref{wdens}
is actually a member of a general family of functionals that have the
form of an LDA, but which are based on a transformed density variable
\begin{align}
  n(\boldsymbol{r}) \rightarrow n(\boldsymbol{r}) f^\text{mLDA}(t(\boldsymbol{r})),
  \label{eq:densmlda}
\end{align}
where $t(\boldsymbol{r})$ is the (dimensionless) reduced kinetic energy density
\begin{equation}
t(\boldsymbol{r})=\frac{\tau(\boldsymbol{r})}{n^{5/3}(\boldsymbol{r})}.\label{eq:t}
\end{equation}
It is easily seen that LDA functionals operating on a density
transformed according to \eqref{densmlda} are exact for the HEG if the
function $f^\text{mLDA}$ reduces to one for the HEG, i.e.
\begin{equation}
  f^\text{mLDA}(C_F) = 1
  \,.
\end{equation}
Because this procedure generates a meta-GGA-type functional without
gradient dependence from a LDA, we will term these functionals
meta-LDAs.

\section{Computational Details\label{sec:Computational-Details}}

The effective density of \eqref{wdens} can be rewritten in the form of
\eqref{densmlda} as
\begin{align}
f(t) & =\left(\frac{t}{C_F}\right)^{3x/5} \,. \label{eq:n-reduced}
\end{align}
The resulting meta-LDA version of the local exchange functional can be
easily rewritten in terms of an enhancement function
\begin{equation}
  F(t;x) = \left[\left(\frac{t}{C_F}\right)^{3x/5}\right]^{4/3}=\left(\frac{t}{C_F}\right)^{4x/5}.\label{eq:Fx-geom}
\end{equation}
The generalization of the Perdew--Wang 1992 correlation
functional\citep{Perdew1992a} is equally trivial; the density used to
evaluate the energy density is merely re-expressed using
\eqref{wdens}. These new functionals have been implemented in version
5.1.0 of the \textsc{Libxc} library of exchange-correlation
functionals.\citep{Lehtola2018} In \textsc{Libxc}, the derivatives of
the functional are evaluated analytically using the \textsc{Maple}
symbolic algebra program, as is the case for all other functionals in
\textsc{Libxc} as well. Combined with a basis set, these derivatives
can be used to minimize the total energy variationally with respect to
the orbital coefficients within a self-consistent field approach; we
refer to \citeref{Lehtola2020} for discussion.

Fully numerical,\citep{Lehtola2019c} fully variational calculations on
closed and partially closed shell atoms from H to Sr were performed
with the finite element method as implemented in the \textsc{HelFEM}
program,\citep{Lehtola2019a} which allows for an efficient approach to
the complete basis set limit.\citep{Lehtola2019e, Lehtola2020c} The
atomic calculations employed five radial elements, yielding 139
numerical radial basis functions which suffice to converge the energy
to better than $\mu E_h$ precision for these systems.

Molecular calculations on the 183 non-multireference molecules in the
W4-17 dataset\citep{Karton2017} were performed with the \textsc{Psi4}
program.\citep{Smith2020} The \textsc{Psi4} calculations employed the
quadruple-$\zeta$ aug-pcseg-3 basis set,\citep{Jensen2001,
  Jensen2002b, Jensen2014} and a $(100,590)$ quadrature grid. Density
fitting\citep{Sambe1975} was used to accelerate the \textsc{Psi4}
calculations; a universal auxiliary basis set was used for this
purpose.\citep{Weigend2008}

\section{Results\label{sec:Results}}

\subsection{Atomic calculations\label{subsec:Atomic-calculations}}

The errors of exchange-only density functional calculations compared
to unrestricted Hartree--Fock (HF) total and exchange energies for
atoms from H to Sr were studied with \textsc{HelFEM}; the reference
unrestricted HF total energies have been recently reported in
\citeref{Lehtola2020c}. Due to the similarity of the results, data is
shown here only for the noble gases Ne, Ar, and Kr in \figref{noble};
the rest of the data can be found in the Supporting Information. In
addition to the self-consistent data, \figref{noble} also shows the
perturbative evaluation of the exchange energy computed on top of the
HF density.

Following \citet{Becke1988a} and \citet{Sun2015} among others, we fit
the parameter $x$ for our meta-LDAs by optimizing the total energy of
the argon atom to the Hartree--Fock reference value, leading to the
choice $x=0.50$. It is noteworthy that in addition to being
quasi-optimal for all systems, $x=0.50$ is also numerically stable for
all the studied atoms. Finally, it also leads to uniformly smaller
errors in the exchange energy than in the LDA and tLDA, which
uniformly underestimate the energy, while LTA grossly overestimates
the energy.

As was already implied above, the self-consistent calculations diverge
for large fractions $x$ of the LTA density. We have analyzed the
instability observed in the calculations; see the Appendix for a
formal analysis. It turns out that the functional form is inherently
unstable for $x>0.625$, since for such values of $x$ the potentials
corresponding to both $n$ and $\tau$ diverge asymptotically to
$-\infty$ for $r \to \infty$. However, it is clear from the results
that the optimal value $x$ for the exchange functional is found at
$x<0.625$.

\begin{figure*}
\subfloat[Ne]{\begin{centering}
\includegraphics[width=0.33\textwidth]{figs/Ne}
\par\end{centering}
}\subfloat[Ar]{\begin{centering}
\includegraphics[width=0.33\textwidth]{figs/Ar}
\par\end{centering}
}\subfloat[Kr]{\begin{centering}
\includegraphics[width=0.33\textwidth]{figs/Kr}
\par\end{centering}
}

\caption{Errors in self-consistent total (blue solid line) and
  exchange (red solid line) energies of Ne, Ar, and Kr, as well as in
  the perturbative exchange energy calculated on top of the HF density
  (dashed red line). The vertical dashed blue line shows the critical
  value $x=0.625$, see main text. The location of the smallest error
  for the self-consistent total and exchange energies are shown as the
  blue and red squares, respectively, and the one for the perturbative
  exchange energy as red diamonds; however, since the optimal value is
  close to $x=1/2$ for all cases, the markers are on top of each
  other.
\label{fig:noble}}
\end{figure*}

\subsection{Molecular calculations\label{subsec:Molecular-calculations}}

\begin{table*}
  \caption{Mean absolute error (MAE) and mean error (ME) in
    atomization energies of the non-multireference part of the
    W4-17 test set, computed in the aug-pcseg-3 basis with density
    fitting and a (100,590) grid. The data is divided into
    exchange-only calculations (\ref{tab:xfunc}), and calculations
    including both exchange and correlation (\ref{tab:xcfunc}).
    See the main text for the legend of the functionals shown. To
    clarify the notation, the used values for $x$ in the meta-LDA
    exchange and correlation functionals are also
    shown.
  }
  \begin{centering}
    \subfloat[Results for exchange-only
      calculations. \label{tab:xfunc}]{\begin{centering}
        \begin{threeparttable}
          \begin{tabular}{llr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}l}
            Functional & $x$ & \multicolumn{2}{c}{MAE (kcal/mol)} & \multicolumn{2}{c}{ME (kcal/mol)}\tabularnewline
            \hline
            LDA exchange &  & 28&966 & -12&015\tabularnewline
            \hline
            hLTA exchange$^{a}$ & $1/2$ & 71&235 & -67&512\tabularnewline
            tLTA exchange & $1/3$ & 47&504 & -35&863\tabularnewline
            qLTA exchange$^{b}$ & $1/4$ & 42&181 & -26&070\tabularnewline
            \hline
            HF &  & 144&848 & -144&848\tabularnewline
            B88 exchange &  & 98&177 & -98&177\tabularnewline
            PBE exchange &  & 87&958 & -87&958\tabularnewline
            \hline
            \hline
          \end{tabular}
          \begin{tablenotes}
          \item [$^{a}$] The data for the exchange-only hLTA
            calculation excludes \ce{CH2NH2} for which the SCF
            procedure did not converge.
          \item [$^{b}$] qLTA is the same as the tLDA of
            \citeauthor{Eich2014}.
          \end{tablenotes}
        \end{threeparttable}
        \par\end{centering}
    } \par\end{centering}
    \begin{centering}
      \subfloat[Results for exchange-correlation
        functionals. \label{tab:xcfunc}]{\begin{centering}
          \begin{threeparttable}
            \begin{tabular}{lllr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}l}
            Functional & $x_{\text{exchange}}$ & $x_{\text{correlation}}$ & \multicolumn{2}{c}{MAE (kcal/mol)} & \multicolumn{2}{c}{ME (kcal/mol)}\tabularnewline
            \hline
            LDA-PW92 &  &  & 79&879 & 79&879\tabularnewline
            \hline
            qLTA-qPW92 & $1/4$ & $1/4$ & 61&089 & 60&897\tabularnewline
            tLTA-tPW92 & $1/3$ & $1/3$ & 50&207 & 49&494\tabularnewline
            hLTA-PW92 & $1/2$ & $0$ & 31&388 & 23&913\tabularnewline
            hLTA-hPW92 & $1/2$ & $1/2$ & 26&907 & 14&088\tabularnewline
            \hline
            B88-P86 &  &  & 19&173 & 18&899\tabularnewline
            PBE-PBE &  &  & 18&028 & 17&052\tabularnewline
            TPSS-TPSS &  &  & 12&427 & 11&180\tabularnewline
            B88-LYP &  &  & 8&176 & 1&714\tabularnewline
            \hline
            \hline
          \end{tabular}
        \end{threeparttable}
          \par\end{centering}
      }
      \par\end{centering}
      \label{tab:Errors-in-atomization}
\end{table*}

The application of the functional to atomization energies
\begin{equation}
E^{\text{at}}=\left[\sum_{\text{atoms }A}E(A)\right]-E(\text{molecule})\label{eq:Eat}
\end{equation}
of the non-multireference part of W4-17 yields the errors
\begin{equation}
\Delta E^{\text{at}}=E^{\text{at}}(\text{DFT})-E^{\text{at}}(\text{W4-17})\label{eq:dEat}
\end{equation}
shown in \tabref{Errors-in-atomization}. Due to the higher cost of the
molecular calculations compared to the atomic ones, the new family of
meta-LDA functionals is only studied at select points which suffice
for the present purposes of showing the proof of concept. The points
at which the meta-LDAs are evaluated are indicated by a prefix to the
name of the exchange and correlation functionals: data are presented
for the LDA exchange functional as qLTA (same as
\citeauthor{Eich2014}'s tLDA), tLTA, and hLTA, which stand for for
$x=1/4$, $x=1/3$, and $x=1/2$, respectively. Data is given both for
exchange-only calculations, and for combinations with the Perdew--Wang
(PW92) correlation functional\citep{Perdew1992a} that also admits
meta-LDA generalizations to qPW92, tPW92 and hPW92 for $x=1/4$,
$x=1/3$, and $x=1/2$, respectively.

For comparison, data is also included for the Perdew--Burke--Ernzerhof
exchange-correlation functional;\citep{Perdew1996, Perdew1997}
combinations of the Becke'88 (B88) exchange
functional,\citep{Becke1988a} with the Perdew'86\citep{Perdew1986,
  Perdew1986b} (P86) and Lee--Yang--Parr\citep{Lee1988} (LYP)
correlation functionals; as well as the
Tao--Perdew--Staroverov--Scuseria (TPSS) exchange-correlation
functional.\citep{Tao2003, Perdew2004}

Starting out with the basics, the table demonstrates the well-known
characteristics of HF and LDA: HF severely underbinds molecules due
to the complete neglect of electronic correlation effects, while LDA
overbinds them. Due to the overbinding, exchange-only LDA calculations
are more accurate than those that explicitly include also correlation
contributions, although LDA exchange by itself is slightly underbinding.
In contrast, while the gradient-corrected exchange functionals yield
bad results if used alone, when they are combined with a good gradient-corrected
correlation functional they achieve good accuracy. Jacob's ladder\citep{Perdew2001}
is also visible in the results: more accurate atomization energies
are obtained in the sequence LDA $\to$ PBE $\to$ TPSS.

Also the meta-LDA functionals interestingly show monotonic behavior.
Going from LDA to qLTA to tLTA and, finally, hLTA in exchange-only
calculations leads to systematically increasing underbinding. The
same effect holds also in the presence of correlation: while LDA-PW92
is greatly overbinding, as was already established above, the overbinding
decreases systematically in the sequence LDA-PW92 $\to$ qLTA-qPW92
$\to$ tLTA-tPW92 $\to$ hLTA-hPW92. Like in the case of the atomic
exchange energies, the half-and-half $x=1/2$ mixture of the electron
density with the $\tau$-based density as in the hLTA-hPW92 functional
yields the best results with a mean absolute error almost three times
smaller than in the original LDA-PW92 calculation. This finding is
underlined by the error histograms shown in \figref{Errors-in-atomization}:
while LDA-PW is consistently overbinding, the errors for hLTA-hPW
are almost symmetric, even though the error scale is still large compared
to established GGA functionals.

\begin{figure*}
\subfloat[LDA-PW]{\begin{centering}
\includegraphics[width=0.5\textwidth]{figs/lda}
\par\end{centering}
}\subfloat[hLTA-hPW]{\begin{centering}
\includegraphics[width=0.5\textwidth]{figs/hlta_hpw}
\par\end{centering}
}

\caption{Error histograms for the atomization energies of the non-multireference
part of W4-17 in the aug-pcseg-3 basis set.\label{fig:Errors-in-atomization}}
\end{figure*}

\section{Summary and Conclusions \label{sec:Conclusions}}

We have proposed a new class of functionals as generalizations of the
established class of local density approximations (LDAs) by including
a fraction $x$ of fictitious density computed from the local kinetic
energy density via a relation derived for the homogeneous electron gas
(HEG). The resulting so-called meta-LDA functionals maintain the
exactness of LDA for the HEG, and are derived from HEG data only
(with the exception of the one parameter $x$ that is fitted to the total
exchange-only energy of the argon atom), but afford much improved
accuracy for inhomogeneous systems, thus forming a new rung on Jacob's
ladder of density functionals in-between LDAs and GGAs. Benchmarks on
both perturbative and self-consistent atomic exchange energies, as
well as molecular atomization energies in the presence of a
correlation functional showed that the half-and-half ratio $x=1/2$
yields quasi-optimal results for both atoms and molecules, almost
fully eliminating the overbinding of LDA and reducing the mean
absolute error in the atomization energies to a third of the original.

Meta-LDAs could also been seen as a better starting point for the
inclusion of an extra dependency in the gradient of the density (as in
a standard GGA), and in the Laplacian of the density and the kinetic
energy density (as in a standard meta-GGA). Due to the extra
flexibility we can expect that these will be better performing than
the parent functionals. For example, the new degree of freedom
introduced with the meta-LDAs could play an important role for, e.g.,
semi-empirical functionals fitted to experimental data. In many of
these cases (see for example on \citerefs{Becke1997, Boese2000,
  Verma2017}), the functionals do not reduce to the LDA for
homogeneous densities, as this would compromise the accuracy of the
functional for other systems. By replacing the standard LDA by a
meta-LDA form in full or in part on could, in principle, obey the
exact condition without compromising accuracy, and at the same time
increase the transferability of the functionals to solids. Of course,
the GGA or meta-GGA enhancement functionals has to be redesigned (or
at least re-optimized) to take the new form into account. Work along
these lines has already started.

\section*{Appendix: instability of the local tau approximation}
The Kohn--Sham electron density is known to behave asymptotically as
$n({\bf r}) \propto \exp(-2\sqrt{-2 \epsilon_\text{HOMO}}r)$ due to
the highest occupied molecular orbital (HOMO) which behaves as
$\psi_\text{HOMO} \propto
\exp(-\sqrt{-2\epsilon_\text{HOMO}}r)$.\cite{Katriel1980} For
simplicity, we will study hydrogenic orbitals of the form
\begin{equation}
  \psi_\sigma ({\bf r}) = \frac {2 \zeta^{3/2} \exp(-\zeta r)}
      {\sqrt{4\pi}} \label{eq:asymptotic}
\end{equation}
to show that the exponentially decaying asymptotic region leads to
problems for $r \to \infty$ for the local tau exchange functional. The
electron density of the wave function in \eqref{asymptotic} is
\begin{equation}
  n_\sigma ({\bf r}) = \frac {\zeta^{3} \exp(- 2\zeta r)}
  {\pi}, \label{eq:den}
\end{equation}
while the kinetic energy density is
\begin{equation}
  \tau_\sigma ({\bf r}) = \frac 1 2 \sum_i |\nabla \psi({\bf r})|^2 =
  \frac {\zeta^5 \exp(-2\zeta r)}{2 \pi}. \label{eq:tau}
\end{equation}
The self-consistent implementation of the meta-LDAs is based on the
potentials $v_n^\sigma$ and $v_\tau^\sigma$, which are defined as the
derivatives of the exchange energy density arising from the
substitution of \eqref{wdens} into \eqref{lda} with respect to
$n_\sigma$ and $\tau_\sigma$, respectively.\cite{Lehtola2020} It is
easy to show using e.g. computer algebra (we used \textsc{Maple 2020}
to obtain these results) that when evaluated on an electron density
and kinetic energy density of the form of \eqref{den, tau}, both of
the potentials $v_n^\sigma$ and $v_\tau^\sigma$ contain a factor of
the form $\exp\left[ (16x-10)r\zeta/15 \right]$, which diverges in the
limit $r\to \infty$ whenever $x > x_\text{crit}$ with the critical
value $x_\text{crit} = 5/8 = 0.625$.

Interestingly, also the choice of a HOMO with a Gaussian form
$\psi_\sigma \propto \exp(-\zeta r^2)$ leads to divergent
potentials---only now of a stronger kind $\exp\left[
  (16x-10)r^2\zeta/15 \right]$---and yields the same critical value
$x_\text{crit}=5/8$. In fact, it can be shown that \emph{all}
asymptotic wave functions of the kind $\psi_\sigma \propto \exp(-\zeta
r^p)$ with $p>0$ lead to divergences of the kind $\exp\left[ (16x-10)
  \zeta r^p /15 \right]$ in $v_n^\sigma$ and $v_\tau^\sigma$. The
total exchange energy, however, is finite in each case.

For $x>x_\text{crit}$ one then has $v_n^\sigma \to -\infty$ and
$v_\tau^\sigma \to -\infty$ for $r \to \infty$, because the potentials
are negative everywhere (as expected for an exchange functional). This
divergence causes convergence problems. Assuming an orthonormal
basis set $\{\chi_\mu\}$, the potentials $v_n^\sigma$ and
$v_\tau^\sigma$ contribute to the Kohn--Sham--Fock matrix
as\cite{Lehtola2020}
\begin{widetext}
\begin{equation}
  F^\sigma_{\mu \nu} \propto \int \left[ v_n^\sigma ({\bf r}) \chi_\mu
    ({\bf r}) \chi_\nu ({\bf r}) + \frac 1 2 v_\tau^\sigma ({\bf r}) \left[
    \nabla \chi_\mu ({\bf r}) \cdot \nabla \chi_\nu ({\bf r}) \right]
    \right] {\rm d}^3r. \label{eq:Ftau}
\end{equation}
\end{widetext}
The tentative physical interpretation of the divergent negative
potentials is that displacing electron density toward $r \to \infty$
would lead to a decrease in the energy. Now, if a Gaussian-type or
Slater-type orbital basis set is employed, $\chi_\mu$ and its gradient
will decay asymptotically as $\exp(-\alpha_\mu r^2)$ or
$\exp(-\zeta_\mu r)$, respectively, where $\alpha_\mu$ and $\zeta_\mu$
are the Gaussian and Slater-type exponent, with analogous expressions
for $\nabla \chi_\nu$. Evaluating \eqref{Ftau} then requires
quadrature of an expression with an exponentially decaying part and an
exponentially increasing part, which is numerically unstable, as the
resulting value may be either small or large. The finite element
calculations with \textsc{HelFEM}, in turn, feature localized basis
functions also at large values of $r$. This leads to exponentially
increasing elements of the Kohn--Sham--Fock matrix, making the
self-consistent field algorithm unstable.

In contrast, the potentials arising in the asymptotic region for
$x<0.625$ decay exponentially (like they do in the local density
approximation), making self-consisistent field calculations stable.

\section*{Acknowledgments}
This work has been supported by the Academy of Finland (Suomen
Akatemia) through project number 311149. Computational resources
provided by CSC -- It Center for Science Ltd (Espoo, Finland) and the
Finnish Grid and Cloud Infrastructure (persistent identifier
urn:nbn:fi:research-infras\-2016072533) are gratefully acknowledged.

\section*{Data Availability}
The data that supports the findings of this study are available within
the article and its supplementary material.

\section*{Supporting Information}

\input{supplement_text.tex}


\section{#1}}
\newcommand{\nsubsection}[1]{\subsection{#1}}
\newcommand{\nsubsubsection}[1]{\subsubsection{#1}}

\else

\newcommand{\ncaption}[1]{\small\caption{#1}}
\newcommand{\nsection}[1]{\vspace{-0.40cm}\section{#1}\vspace{-0.3cm}}
\newcommand{\nsubsection}[1]{\vspace{-0.4cm}\subsection{#1}\vspace{-0.15cm}}
\newcommand{\nsubsubsection}[1]{\vspace{-0.3cm}\subsubsection{#1}\vspace{-0.15cm}}

\fi

\input{math_commands.tex}

\hypersetup{
 colorlinks,
  linkcolor={green!80!black},
  citecolor={red!70!black},
  urlcolor={blue!70!black}
}

\begin{document}
\date{}

\title{\Large \bf Towards Robust LiDAR-based Perception in Autonomous Driving: General Black-box Adversarial Sensor Attack and Countermeasures}

\author{
{\rm Jiachen Sun}\\
University of Michigan\\
\href{mailto:jiachens@umich.edu}{jiachens@umich.edu}
\and
{\rm Yulong Cao}\\
University of Michigan\\
\href{mailto:yulongc@umich.edu}{yulongc@umich.edu}
\and
{\rm Qi Alfred Chen} \\
UC Irvine\\
\href{mailto:alfchen@uci.edu}{alfchen@uci.edu}
\and
{\rm Z. Morley Mao}\\
University of Michigan\\
\href{mailto:zmao@umich.edu}{zmao@umich.edu}
}

\maketitle

\input{sections/Abstract.tex}
\input{sections/Introduction.tex}
\input{sections/Background.tex}
\input{sections/Threat_Model.tex}
\input{sections/Limitations.tex}
\input{sections/Vulnerability.tex}
\input{sections/Attack.tex}
\input{sections/Ad_Defense.tex}
\input{sections/Ml_Defense.tex}
\input{sections/Discussion.tex}
\input{sections/Related_Works.tex}
\input{sections/Conclusion.tex}

\renewcommand{\bibfont}{\footnotesize}
\bibliographystyle{abbrv}

\subsubsection{Hierarchies}

\nsubsubsection{Hierarchy Design}

To achieve both robustness and efficiency, CARLO hierarchically integrates FSD and LPD. In the first stage, CARLO accepts the detected bounding boxes and leverages LPD to filter the unquestionably fake and valid vehicles by two thresholds (\S\ref{ltd}). The rest bounding boxes are uncertain and will be further fed into FSD for final checking. CARLO achieves around 8.5 ms mean processing time for each vehicle. The entire algorithm of CARLO is detailed in Appendix \ref{ap:algorithm}.

\nsubsection{CARLO Evaluation}

\textbf{Experimental setup.} We evaluate the defense performance of CARLO on the KITTI trainval and test sets. We apply all the attack traces from $\mathcal{K},\mathcal{R}$ to all point cloud samples at target locations (5-8 meters in front of the victim), and feed them into three CARLO-guarded models $\mathsf{CARLO}(\mathcal{M}(\cdot))$. We also evaluate CARLO against Adv-LiDAR~\cite{cao2019adversarial} on Apollo 5.0. The defense goal is to successfully detect the spoofed fake vehicles from the output bounding boxes without hurting the original performance (\textit{i.e.,}\xspace AP) of the target models.

\begin{figure}[!tp]
\vspace{-0.4cm}
\centering
\subfigure[CARLO-guarded Apollo 5.0.]{
\begin{minipage}[t]{0.49\linewidth}
\centering
\includegraphics[width=4cm]{images/usenix_ad_apollo.pdf}
\end{minipage}
}%
\subfigure[CARLO-guarded PointPillars.]{
\begin{minipage}[t]{0.49\linewidth}
\centering
\includegraphics[width=4cm]{images/usenix_ad_pointpillars.pdf}
\end{minipage}%
}%
\vskip 0pt
\vspace{-0.3cm}
\subfigure[CARLO-guarded PointRCNN.]{
\begin{minipage}[t]{0.49\linewidth}
\centering
\includegraphics[width=4cm]{images/usenix_ad_pointrcnn.pdf}
\end{minipage}
}%
\subfigure[Precision and recall of CARLO.]{
\begin{minipage}[t]{0.49\linewidth}
\centering
\includegraphics[width=4cm,height=2.4cm]{images/usenix_ad_pr.png}
\end{minipage}
}%

\centering
\vspace{-0.2cm}
\caption{Attack success rates ($ASR$s) of proposed black-box spoofing attacks on three CARLO-guarded models.}
\label{fig:asr_rhd}
\vspace{-0.5cm}
\end{figure}

\textbf{Evaluation metrics.} We evaluate the performance of CARLO in two aspects, which are the $ASR$ of the CARLO-guarded models, and the precision and recall of CARLO itself. $ASR$ directly relates to the defense performance, while the precision and recall of CARLO reflect whether it will harm the original AP of target models. We test the $ASR$ on the validation set and test set since the distributions are estimated from the training set. We only evaluate the precision and recall of CARLO on the validation set as we do not have the ground-truth for the test set.

Figure \ref{fig:asr_rhd} (a-c) shows the $ASR$ of three CARLO-guarded models. As shown, CARLO reduces the $ASR$ from more than 95\% to below 9.5\% with the maximum attack capability, and reduce the mean $ASR$ to around 5.5\%. \CR{We observe that the remaining 5.5\% comes from the detection errors (\textit{i.e.,}\xspace the detected bounding box of the fake vehicle cannot match well with the ground-truth) that shift the $f'$ and $g'$ to the distribution of valid vehicles.} The errors occur randomly in the point clouds so that it is hard for adversaries to utilize. The recall in Figure \ref{fig:asr_rhd} (d) reaches around 95\% in all targets models which also validate the results in Figure \ref{fig:asr_rhd} (a-c).

\begin{table}[t]
    \captionsetup{font=small}
    \footnotesize
    \caption{PointPillars' and PointRCNN's APs (\%) of 3D car detection on the KITTI validation set. ``Mod.'' refers to the Moderate category introduced in \S\ref{background_kitti}; ``Original'' refers to the original performance of two models; ``Attack'' refers to the performance after spoofing attacks; ``CARLO'' refers to the performance after CARLO applied.}
    \vspace{-0.1in}
    \centering
    \begin{tabular}{|c|c|c|c|c|c|c|}
    \hline
    \multirow{2}*{Model} &  \multicolumn{3}{ c |}{PointPillars}  &  \multicolumn{3}{ c |}{PointRCNN} \\
    \cline{2-7}
     &  Easy & Mod. & Hard  &  Easy & Mod. & Hard \\
    \hline
    \hline
    Original & 86.56 & 76.87 & 72.09 &  88.80 & 78.58 & 77.64 \\
    \hline
    Attack &  74.06 & 56.69 & 53.98 & 84.51 & 71.17 & 68.06 \\
    \hline
    CARLO &  86.57 & 78.60 & 73.55 & 88.91 & 78.61 & 77.63 \\
    \hline
    \end{tabular}
    \label{tb:model_og_performance}
    \vspace{-0.5cm}
\end{table}

Besides delivering satisfactory defense performance, CARLO barely introduces misdetections (\textit{i.e.,}\xspace false negatives) to the models. Figure \ref{fig:asr_rhd} (d) shows that the precision of anomaly detection reaches at least 99.5\% for all target models. We manually verify the 0.5\% misdetections, and find they are all vehicles at least 40 meters away from the AV, which will not affect its immediate driving behavior. Table \ref{tb:model_og_performance} also shows that the AP will slightly increase after CARLO being applied to the original model because the original model has internal false positives such as detecting a flower bed as a vehicle. CARLO detects some of those false positives generated by the original model.

\nsubsubsection{Defense against White-box and Adaptive Attacks}

To further evaluate CARLO against white-box attacks, we first leverage Adv-LiDAR to generate adversarial examples that fool Apollo 5.0 and test whether they can succeed in attacking CARLO-guarded Apollo 5.0. Figure \ref{fig:white_box_rhd} demonstrates that CARLO can effectively defend Adv-LiDAR, where the $ASR$ drops from more than 95\% to below 5\% consistently. We observe that the defense effects are better than the results shown in Figure~\ref{fig:asr_rhd} (a). We find out that Adv-LiDAR tends to translate the attack traces to a slightly higher place along $z$ axis. Such translations will isolate the adversarial examples in the point cloud, making themselves easier to be detected by CARLO.

We also try our best efforts to evaluate CARLO on the adaptive attacks. We assume attackers are aware of the CARLO pipelines and utilize Adv-LiDAR to break CARLO's defense. Due to the sensor attack capability, adversaries have limited ability to modify the absolute free space ($\sum_{\,c \in \mathsf{B}}\mathds{1} \cdot FS(c)$) in Equation \ref{eq:fsd}. However, attackers can try to shrink the volume of the bounding box ($|\mathsf{B}|$) to shift the distribution of $f'$, since it is controlled purely by models. Therefore, the attack goal is to spoof a vehicle at target locations at the same time, minimize the size of the bounding box. We formulate the loss function and follow Adv-LiDAR~\cite{cao2019adversarial} to utilize a global transformation matrix $H(\theta,\tau)$ (Equation~\ref{eq:trans}) for solving the optimization problem:
\CR{
\begin{equation}
\begin{split}

    & \min_{\theta,\tau} \qquad {\mathcal{L}} (x \oplus  V \cdot H(\theta,\tau)^{T}) + \lambda \cdot {\mathcal{V}}_{B} ( V \cdot H(\theta,\tau)^{T})

\end{split}
\label{eq:adaptive_1}
\end{equation}
where ${\mathcal{L}}(\cdot)$ is the loss function defined in Adv-LiDAR~\cite{cao2019adversarial}, ${\mathcal{V}}_{B}(\cdot)$ is the volume of the target bounding box $B$, and $\lambda$ is a hyper-parameter.} Figure~\ref{fig:adaptive_rhd} shows that such adaptive attacks cannot break CARLO, either. We attribute the reason to $H(\theta,\tau)$ that holistically modifies the spoofed points so that it can barely change the size of the bounding box.

\begin{figure}[t]
  \vspace{-0.3cm}
  \begin{minipage}[t]{0.49\linewidth}
    \centering
    \includegraphics[width=\linewidth]{images/usenix_whitebox_ad.pdf}
    \caption{Attack success rates ($ASR$s) of Adv-LiDAR on Apollo 5.0 and CARLO-guarded model.}
    \label{fig:white_box_rhd}
  \end{minipage
  \hfill
  \begin{minipage}[t]{0.48\linewidth}
    \centering
    \includegraphics[width=\linewidth]{images/usenix_ad_adaptive.pdf}
    \caption{Attack success rate ($ASR$) of the adaptive attack on CARLO-guarded Apollo 5.0.}
    \label{fig:adaptive_rhd}
  \end{minipage}
  \vspace{-0.5cm}
\end{figure}

\subsection{Digital Spoofing Attack}

\emph{1. Constructing original attack traces.} As demonstrated in \S\ref{vulnerability_id}, occluded or distant vehicles' point sets that meet the sensor attack capability can be utilized to spoof front-near vehicles. Therefore, our methodology attempts to closely represent realistic physical attacks using traces from real-world datasets (\textit{e.g.,}\xspace KITTI). In order to test different sensor attack capability, we extract occluded vehicles' point sets with varying numbers of points (5-200 points) from the KITTI validation set. Furthermore, we take 10 points as interval, and divide the extracted point sets into 20 groups per their number of points (The first group contains traces with the number of points from 0 to 10, and the second group contains traces with the number of points from 10 to 20, \textit{etc.}\xspace). We then randomly pick five traces in each group forming a small dataset $\mathcal{K}$ containing 100 point sets.

Besides collecting existing real-world traces, the identified vulnerability also supports adversaries in generating customized attack traces, which are more efficient for pipelining the attack process. We leverage ray-casting techniques to generate customized attack traces. More specifically, we utilize a 3D car mesh and implement a renderer~\cite{cao2019adversarial2} simulating the function of a LiDAR sensor that probes the car mesh by casting lasers. By doing so, we can render the car mesh's point cloud. We further simulate different occlusion and sparsity patterns on the car mesh to fit the sensor attack capability, as shown in Figure \ref{fig:render}. Similar to $\mathcal{K}$, we collect rendered point clouds with different numbers of points by using different postures and occlusion patterns. We also follow the same procedure to build a small dataset $\mathcal{R}$ containing 100 rendered point sets. More figures of $\mathcal{R}$ are shown in Appendix \ref{ap:figure}.

\emph{2. Spoofing original attack traces at target locations.} To trigger severe security and safety consequences, adversaries need to inject the constructed attack traces at target locations in the point cloud.  We consider spoofing $\mathcal{K}$ and $\mathcal{R}$ in both digital and physical environments. For digital spoofing, we make sure the injection of attack traces meets the sensor attack capability $\mathcal{A}$ and real-world requirements. We follow the high-level formulation in Adv-LiDAR~\cite{cao2019adversarial} utilizing a global transformation matrix $H(\theta,\tau)$ (Equation \ref{eq:trans}) to translate the attack traces (\textit{i.e.,}\xspace ${V'}^{T} = H(\theta,\tau)\cdot V^{T}$, where $V \in \mathcal{K} \cup \mathcal{R}$). Here the translation interprets the attack capability ($\mathcal{A}$) in terms of modifying the \emph{azimuth} and \emph{distance} of attack traces. We further calibrate each point in the translated attack trace to its nearest laser ray's direction and prune the translated attack trace to fit the attack capability (\textit{i.e.,}\xspace $V' \in \mathcal{A}$). Finally, we merge the attack trace with the pristine point cloud according to the physics of LiDAR. We feed the modified point cloud samples containing the attack traces into three target models. For physical spoofing, we program attack traces from $\mathcal{R}$ as input to the function generator so that we can control the spoofed points and launch the spoofing attack~\cite{shin2017illusion} in our lab. We further collect the physical attack traces and feed them into target models. Due to the limitation of our attack devices, we only conduct preliminary physical spoofing experiments. More details of physical spoofing can be found in \S\ref{practi}. It is worth noting that such limitations do not hurt the validity of our attack model (\S\ref{threat}) since the attack capability $\mathcal{A}$ is adopted from Adv-LiDAR~\cite{cao2019adversarial}, in which has been demonstrated in the real world.

\nsubsection{Attack Evaluation and Analysis}
\label{attack_analysis}
We perform large-scale evaluations on our proposed black-box attack in terms of effectiveness and robustness.

\textbf{Experimental setup.} \CR{The evaluations are performed on the KITTI trainval and test sets (introduced in~\S\ref{background_kitti}), which are collected in the physical world. As mentioned before, limited by our attack devices, we leverage $\mathcal{K},\mathcal{R}$ to launch \textit{digital spoofing} attacks. We also utilize attack traces ($\mathcal{S}$) generated by the sensor-level spoofing attack (\S\ref{limitations_existing}) as a baseline. $\mathcal{S}$ is collected from blindly \textit{physical spoofing} attacks on a real Velodyne VLP-16 PUCK LiDAR~\cite{vpl16}. We further inject all the attack traces from above three constructed datasets into the KITTI point clouds at front-near locations (\textit{i.e.,}\xspace 5-8 meters in front of the victim AV) to test their effectiveness.}

\textbf{Evaluation metrics.} Object detection models often have default thresholds for confidence scores to filter out detected objects with low confidence (potential false positives). We leverage the default thresholds used by three target models to measure the attack success rate ($ASR$). We label an attack successful as long as the model detects a vehicle at the target location whose confidence score exceeds the default threshold:
\begin{equation}
    ASR = \frac{\mathrm{\#\; of\; successful\; attacks}}{\mathrm{\#\; of\; total\; point\; cloud\; samples}}
    \label{eq:asr}
\end{equation}

Besides the default threshold, we also define a new metric that leverages multiple thresholds to evaluate LiDAR spoofing attacks. The corresponding definitions and evaluations are described in Appendix \ref{ap:attack_eval}, which provide insights that point-wise features appear to be more robust than voxel-based features.

\nsubsubsection{Attack Effectiveness}
\label{attack_effe}

Figure \ref{fig:asr_all} shows the $ASR$ of the digital spoofing attack with different attack capabilities (\textit{i.e.,}\xspace number of points). As expected, the $ASR$ increases with more spoofed points. The $ASR$s are able to universally achieve higher than 80\% in all target models with more than 60 points spoofed, and it also stabilizes to around 85\% with more than 80 points spoofed. Notably, the attack traces from $\mathcal{R}$ can achieve comparable $ASR$ with $\mathcal{K}$ on all target models, which demonstrates that adversaries can efficiently leverage a customized renderer to generate attack traces (Figure \ref{fig:render}). Such rendered traces can be directly programmed into hardware for physical spoofing attacks (Appendix \ref{ap:Spoof}). Interestingly, $\mathcal{S}$ achieves much higher $ASR$ on Apollo 5.0, indicating that BEV-based features are less robust to spoofing attacks than the other two categories, which could be attributed to the information loss of feature encoding from BEV.

\nsubsubsection{Robustness Analysis}

We analyze the robustness of the proposed attack to variations of attack traces $V'$ and the average precision (AP) of target models $\mathcal{M}$. We also evaluate the attack robustness against state-of-the-art defense strategies~\cite{xu2017feature,yang2019menet} designed for image-based adversarial attacks. We find that spoofed traces with around 60 points to trigger major changes in $ASR$. Note that Cao \textit{et al.}\xspace also utilized spoofed traces with 60 points for analysis~\cite{cao2019adversarial}. Therefore, we use attack traces with $(60,70]$ points from $\mathcal{R}$ for the robustness analysis.

\textbf{Robustness to variations in attack traces}. First, we apply a scaling matrix $S$ to the attack traces $V'$ with different-level randomness to simulate the inaccuracy of sensor attack:
\begin{equation}
\begin{split}
& \quad {V''}_{\wi} = {V'}_{\wi} \\
\begin{bmatrix}
  {V''}_{\mathbf{w_x}}\\
  {V''}_{\mathbf{w_y}}\\
  {V''}_{\wz}\\
\end{bmatrix} & = \begin{bmatrix}
  s & 0 & 0\\
  0 & s & 0\\
  0 & 0 & s\\
\end{bmatrix}
\cdot
\begin{bmatrix}
  {V'}_{\mathbf{w_x}}\\
  {V'}_{\mathbf{w_y}}\\
  {V'}_{\wz}\\
\end{bmatrix}
\end{split}
\end{equation}
where $s$ subjects to a uniform distribution $U(1-\epsilon,1+\epsilon)$.
We use the mean $l$-$2$ norm to measure the distance between $V''$ and $V'$. Figure \ref{fig:robust_1} shows the $ASR$ drops with larger $l$-$2$ distances which is expected. However, as shown, the $ASR$ still reaches around 70\% while the distance is around 10~cm. We also observe that the $ASR$ for PointRCNN drops faster than for PointPillars and Apollo 5.0, which also validates that point-wise features are arguably more robust than voxel-based and BEV-based features (detailed in Appendix \ref{ap:attack_eval}).

\begin{figure}[!t]
  \begin{minipage}[t]{0.478\linewidth}
    \centering
    \includegraphics[width=\linewidth]{images/usenix_robustness1.pdf}
    \vspace{-0.5cm}
    \caption{Attack robustness to variations in generated attack traces from $\mathcal{R}$.}
    \label{fig:robust_1}
  \end{minipage
  \hfill
  \begin{minipage}[t]{0.49\linewidth}
    \centering
    \includegraphics[width=\linewidth]{images/usenix_robustness2.pdf}
    \vspace{-0.5cm}
    \caption{Attack robustness to variations in target models' performance.}
    \label{fig:robust_2}
  \end{minipage}
  \vspace{-0.3cm}
\end{figure}

\textbf{Robustness to variations in model performance}. To understand the relationship between $ASR$ and the original performance of models (\textit{i.e.,}\xspace AP), we first extract the intermediate models when we trained PointPillars and PointRCNN. We then try to launch attacks on these models. Surprisingly, we find that the $ASR$ increases with higher AP\footnote{\CR{We evaluate the $ASR$s
until the training procedures (\textit{i.e.,}\xspace APs) converge on both models.}} (Figure \ref{fig:robust_2}), which implies that a model with better performance could be more vulnerable to such attacks. Our results indirectly demonstrate that the identified vulnerability could be attributed to an ignored dimension (\textit{i.e.,}\xspace occlusion patterns) by current models. Since the models do not notice such a hidden dimension, they will be overfitted to be more vulnerable during training.

\CR{\textbf{Robustness to adversarial training}. Adversarial training is not rigorously applicable because it targets classification models, and requires norm-bounded perturbations to make the optimization problem tractable~\cite{madry2017towards}. In contrast, our study targets 3D object detection models, and the proposed attack is constrained by the sensor attack capability ($\mathcal{A}$), which does not fit any existing norm-bounded formulations. Thus, we perform this robust analysis in an empirical setting. Specifically, we generate another 100 attack traces with 60 points using the customized renderer and randomly inject two of them into each point cloud sample in the KITTI training set at areas without occlusions. We further train PointPillars and PointRCNN on this modified dataset and evaluate our proposed attack using the same 60-point attack traces with \S\ref{attack_effe} on them. We observe that the $ASR$s drop from 83.6\% to 70.1\% and 88.3\% to 79.7\% on PointPillars and PointRCNN, respectively, on the KITTI validation set. However, the ``Hard'' category's original detection performance has significant degradation of over 10\% on both models. Our results empirically show that current LiDAR-based perception model designs cannot learn the occlusion information correctly. The slight drop of the $ASR$s comes from the under-fitting effect of existing occluded vehicles (\textit{i.e.,}\xspace significant AP degradation), which is not acceptable in real AD systems.}

\textbf{Robustness to randomization-based defenses}. We leverage state-of-the-art image-based defenses: feature squeezing~\cite{xu2017feature} and ME-Net~\cite{yang2019menet} to test the attack robustness on Apollo 5.0 since it has similar pipelines with image-based models. We demonstrate that none of them can defend the black-box spoofing attack without hurting the original AP. More details can be found in Appendix \ref{ap:attack_eval}.

\subsubsection{BEV-based 3D Object Detection}

\textbf{1. Bird's-eye view (BEV)-based 3D object detection.} Due to the remarkable progress made in 2D image recognition tasks, a large number of existing works~\cite{apollo,yang2018pixor,liang2018deep,meyer2019lasernet} attempt to transform LiDAR point clouds into the 2D structure for 3D object detection in AD systems. Most state-of-the-art methods~\cite{apollo,yang2018pixor,liang2018deep} conduct the transformation by projecting point clouds into the top-down view, also known as the BEV, and utilize
convolutional neural networks (CNNs) to perform the final detection. Figure \ref{fig:sota_models} (a) shows the architecture of \emph{Apollo 5.0}\footnote{In this paper, we use ``Apollo 5.0'' to denote the Baidu Apollo 5.0 model.}, \CR{an industry-level BEV-based model,} that has six hard-coded feature maps in the BEV and follows a UNet-like~\cite{ronneberger2015u} pipeline to output the grid-level confidence score. The final stage heuristically clusters the grids that belong to the same object.

\begin{figure}[t!]
\centering
  \vspace{-0.2cm}
  \includegraphics[width=\linewidth]{images/usenix_existing_model_new.pdf}
  \vspace{-0.3cm}
  \vspace{-0.15in}
  \caption{State-of-the-art LiDAR-based perception models.}
  \label{fig:sota_models}
  \vspace{-0.5cm}
\end{figure}

\textbf{2. Voxel-based 3D object detection.} VoxelNet~\cite{zhou2018voxelnet} is the first model that slices the point clouds into voxels and extracts learnable features by applying a PointNet~\cite{qi2017pointnet} to each voxel, after which a 2D convolutional detection layer is applied in the final stage. Many recent works~\cite{lang2019pointpillars,yan2018second,wang2019voxel,lehner2019patch} adopt this voxel-based architecture and achieve state-of-the-art performance~\cite{kitti_3d}. Figure \ref{fig:sota_models} (b) shows the architecture of \emph{PointPillars} that creatively voxelizes the point cloud into pillars (a representation of point clouds organized in vertical columns) to enhance the efficiency and follows the general design of voxel-based detection architectures. \CR{Notably, PointPillars is adopted by Autoware~\cite{autoware}, an industry-level AV platform.}

\textbf{3. Point-wise 3D object detection.} Instead of transforming point clouds to regular 2D structures or voxels for feature extraction, recent studies propose to directly operate on point clouds for 3D object detection~\cite{shi2019pointrcnn,yang2019std,chen2019fast,shi2019part} and achieve the state-of-the-art performance. Most existing works in this category use a classic two-stage architecture similar to Faster RCNN~\cite{ren2015faster} in 2D object detection. The first stage is responsible for generating high-quality region proposals in the 3D space. Based on these proposals, the second stage regresses the bounding box parameters and classifies the detected objects. As shown in Figure \ref{fig:sota_models} (c), \emph{PointRCNN} adopts a bottom-up method that generates point-wise region proposals in the first stage and regresses these proposals in the later stage.

\nsubsubsection{KITTI Benchmark}
\label{background_kitti}

KITTI~\cite{Geiger2013IJRR} is a popular dataset for benchmarking AD research, of which the point cloud data are by design divided into a \textbf{trainval set} containing 7481 samples and a \textbf{test set} containing 7518 samples. We follow the methodology by Chen \textit{et al.}\xspace to split the trainval set to a \textbf{training set} (3712 samples) and a \textbf{validation set} (3769 samples) for better experimental studies~\cite{chen20153d}. KITTI evaluates 3D object detection performance by average precision (AP) using the PASCAL~\cite{Everingham10} criteria and requires a 3D bounding box overlap (IoU) over 70\% for car detection. KITTI also defines objects into three difficulty classes: Easy, Moderate, and Hard~\cite{kitti_3d}. The difficulties correspond to different occlusion and truncation levels. We train PointPillars and PointRCNN on the training set, and Table \ref{tb:model_og_performance} shows their APs evaluated on the validation set. We utilize the publicly released Apollo 5.0 model since it has its own labeling, which is incompatible with KITTI. \CR{In this work, we target car detection on the KITTI benchmark as the APs of pedestrian and cyclist detection are not yet satisfactory. However, our methodology can be generalized to other categories.}

\nsubsection{LiDAR Sensor and Spoofing Attacks}
\label{background_spoof}

\textbf{LiDAR sensor}. A LiDAR instrument measures the distance to surroundings by firing rapid laser pulses and obtaining the reflected light with a sensor. Since the speed of light is constant, the accurate distance measurements can be derived from the time difference between laser fires and returns. By firing laser pulses at many predetermined vertical and horizontal angles, a LiDAR generates a point cloud that can be used to make digital 3D representations of surroundings. Each point in a point cloud contains its $xyz$-$i$ information, corresponding to its location and the intensity of the captured laser return.

\nsubsubsection{Sensor-level LiDAR Spoofing Attack}
\label{background_direct_spoof}

In the context of sensors, a spoofing attack is the injection of a deceiving physical signal into a victim sensor~\cite{198480}. Since they share the same physical channels, the victim sensor accepts the malicious signal, trusting it as legitimate. Prior works~\cite{petit2015remote,shin2017illusion} have shown that LiDAR is vulnerable to laser spoofing attacks. Specifically, Petit~\textit{et al.}\xspace showed the feasibility to relay LiDAR laser pulses from other locations to inject fake points into the point cloud~\cite{petit2015remote}. Shin~\textit{et al.}\xspace further improved the attack to control fake points at different locations in the point cloud, even very close to the victim vehicles~\cite{shin2017illusion}.

\nsubsubsection{Adv-LiDAR: Model-level LiDAR Spoofing Attack}
\label{background_adv_spoof}
Besides directly spoofing fake points into LiDAR point clouds, a recent study proposes Adv-LiDAR that uses adversarial machine learning to not only spoof a set of fake points into the point cloud but also manage to deceive the LiDAR-based perception model~\cite{cao2019adversarial}. The authors formulate the attack on Apollo 2.5\footnote{Apollo 2.5 was the latest version when Adv-LiDAR~\cite{cao2019adversarial} was published. In this work, we target Apollo 5.0, the latest version at the time of writing.} as an optimization problem:
\begin{equation}
\begin{split}

    & \min \qquad {\mathcal{L}} (x \oplus t';\mathcal{M}) \qquad

    \\ & \text{s.t.} \qquad t' \in  \{ \Phi(T') \,|\, T' \in \mathcal{A} \} \And x = \Phi(X)
\end{split}
\label{eq:adv_lidar}
\end{equation}
where $X$ is the pristine point cloud and $x$ represents the hard-coded feature maps in Apollo (\S\ref{background_percpetion}). $\Phi(\cdot)$ is the pre-processing function for crafting the feature maps. $T'$ and $t'$ are the spoofed point cloud and its corresponding feature maps, respectively. $\mathcal{A}$ stands for the sensor attack capability, and $\oplus$ merges the pristine and adversarial feature maps.

The attack goal is to spoof a fake vehicle right in front of the victim AV that leads to safety issues, \CR{and the success condition is that the confidence score of the optimized spoofed points ($T'$) exceeds the default threshold so that Apollo 2.5 ($\mathcal{M}$) will detect $T'$ as a valid vehicle.} The authors formulate the sensor attack capability ($\mathcal{A}$) for general LiDAR spoofing attacks and design a specific loss function (${\mathcal{L}}$) and a merging function ($\oplus$) for Apollo 2.5 ($\mathcal{M}$). By strategically controlling the spoofed points, Adv-LiDAR achieves around 75\% attack success rate towards Apollo 2.5 and is considered as the state-of-the-art LiDAR spoofing attack.

\subsubsection{Existing View Fusion Frameworks\yl{Does not make sense to have only one subsection.}}

Besides the models that only take FV as input, several studies~\cite{chen2017multi,zhou2019end,ku2018joint}  present fusion-based architectures for LiDAR-based perception that utilize the combinations of data from different sensors and views as input. MV3D~\cite{chen2017multi} is a classic fusion-based design that takes both LiDAR point clouds and RGB images as input and predicts 3D bounding boxes, where the point cloud is projected into multi-views (\textit{i.e.,}\xspace FV and BEV) for feature encoding. Zhou \textit{et al.}\xspace recently proposed multi-view fusion (MVF), which combines FV with 3D representations~\cite{zhou2019end}. MVF builds on top of PointPillars. Instead of only voxelizing points in 3D, MVF also voxelizes the point cloud into FV frustums and integrates the two voxels' features based on coordination mapping in the 3D space.

To better understand the robustness of fusion-based architectures, we reproduce MV3D and MVF based on PointPillars. For MV3D, we ignore the RGB images, and take the FV and BEV as the model input since we focus on LiDAR-based perception. We use a VGG-16~\cite{Simonyan15} for FV feature learning in MV3D. Figure \ref{fig:existing_fusion} shows the architectures we adopt and reproduce. We train the two reproduced models on the KITTI training set and evaluate them on the KITTI validation set. As Table~\ref{tb:fusion_performance} shows, the FV-augmented models can achieve comparable performance than the original PointPillars. The reproduced results also align well with the evaluations in~\cite{zhou2019end,chen2017multi}.

\begin{table}[thp!]
\captionsetup{font=small}
\vspace{-0.25cm}
    \footnotesize
    \caption{MV3D-PointPillars' and MVF-PointPillars' APs (\%) of 3D car detection on the KITTI validation set.}
    \vspace{-0.1in}
    \centering
    \begin{tabular}{|c|c|c|c|}
    \hline
    \multirow{2}*{Model} &  \multicolumn{3}{ c |}{Car Detection} \\
    \cline{2-4}
     &  Easy & Moderate & Hard \\
    \hline
    \hline
    MV3D-PointPillars & 85.67 & 77.12 & 71.65 \\
    \hline
    MVF-PointPillars &  86.77 & 79.15 & 75.72 \\
    \hline
    \end{tabular}
    \label{tb:fusion_performance}
    \vspace{-0.3cm}
\end{table}

We then evaluate their robustness against our proposed black-box attack. The experimental setups are identical to the settings in \S\ref{attack_analysis}. Figure \ref{fig:asr_fusion} shows that the $ASR$ of reproduced models are as high as the original PointPillars. It indicates that existing view fusion-based architectures both cannot help with defending against LiDAR spoofing attacks, although they provide marginally gain on the AP. We further perform ablation studies and find that the BEV (or 3D) features dominate the model decisions (elaborated in Appendix~\ref{ap:ablation}). Since the identified vulnerability exists in 3D space, the models are still vulnerable to LiDAR spoofing attacks.

\begin{figure}[tp]
\vspace{-0.3cm}
\centering

\subfigure[$ASR$ of the reproduced MV3D.]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=0.95\linewidth]{images/usenix_mv3d_asr.pdf}
\end{minipage}%
}%
\subfigure[$ASR$ of the reproduced MVF.]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=0.95\linewidth]{images/usenix_mvf_asr.pdf}
\end{minipage}%
}%
\centering
\vspace{-0.2cm}
\vspace{-0.1in}
\caption{Attack success rates ($ASR$s) of proposed black-box spoofing attacks on MV3D and MVF models.}
\label{fig:asr_fusion}
\vspace{-0.6cm}
\end{figure}

\nsubsection{Sequential View Fusion}
\label{svf}

The insights drawn from existing view fusion schemes show that existing fusion designs cannot provide better robustness compared to the original models. The 3D (or BEV) representation dominates the model leaving the FV representation not critical in the end-to-end architectures.

Based on the above understandings, we propose a new view fusion schema called sequential view fusion (SVF). SVF comprises of three modules (Figure \ref{fig:svf}), which are: 1) semantic segmentation: a semantic segmentation network that utilizes the FV representation to computes the point-wise confidence scores (\textit{i.e.,}\xspace the probability that one point belongs to a vehicle). 2) view fusion: the 3D representation is augmented with semantic segmentation scores. 3) 3D object detection: a LiDAR-based object detection network that takes the augmented point clouds to predict bounding boxes. Instead of leaving the models to learn the importance of different representations by themselves, we attach a semantic segmentation network to the raw FV data. By doing so, we enforce the end-to-end learning to appreciate the FV features, so that the trained model will be resilient to LiDAR spoofing attacks.

\textbf{Semantic segmentation.} The semantic segmentation networks accept the FV represented point clouds and associate each point in FV with a probability score that it belongs to a vehicle. These scores provide aggregated information on the FV representation. Semantic segmentation over FV has several strengths. First, as mentioned before, the FV representation is noisy because of the nature of LiDAR. Compared to 3D object detection or instance segmentation, which is intractable over FV, semantic segmentation is an easier task as it does not need to estimate object-level output. Second, there are extensive studies on semantic segmentation over FV represented point clouds~\cite{wu2018squeezeseg,wang2018pointseg,biasutti2019lu}, and the segmentation networks achieve much more satisfactory results than the 3D object detection task over FV.

In our implementation, we adopt the high-level design in LU-Net~\cite{biasutti2019lu}. It is worth noting that the end-to-end SVF architecture is agnostic to the semantic segmentation module.

\textbf{View fusion.} The fusion module re-architects existing symmetric designs which integrate the 3D representation with the confidence scores generated by the semantic segmentation module. Specifically, we use Equation \ref{eq:mapping} for mapping between $\vec{p}=(x,y,z)$ and $\vec{p}_{FV}(r,c)$, and augment each $\vec{p}$ with the point-wise confidence score from its corresponding $\vec{p}_{FV}$.

\textbf{3D object detection.} SVF is also agnostic to the 3D object detection module. In this paper, we utilize PointPillars and PointRCNN in our implementation. Most of the models introduced in \S\ref{background} can fit into the end-to-end SVF architecture.

\nsubsection{SVF Evaluation}
\label{svf_eval}

\textbf{Experimental setup.} We train SVF-PointPillars and SVF-PointRCNN on the KITTI training set. The setup of robustness analysis against LiDAR spoofing attacks is identical to the settings in \S\ref{attack_analysis}. We also try to evaluate SVF against Adv-LiDAR~\cite{cao2019adversarial} on Apollo 5.0 and the adaptive attacks.

\begin{figure}[tp]
\vspace{-0.2cm}
\centering
  \includegraphics[width=\linewidth]{images/SVF_models.pdf}
  \vspace{-0.25in}
  \caption{Sequential view fusion (SVF) architecture.}
  \label{fig:svf}
  \vspace{-0.5cm}
\end{figure}

\textbf{Evaluation metrics.} We evaluate the AP of SVF-PointPillars and SVF-PointRCNN on the KITTI validation set, and leverage $ASR$ to test their robustness against LiDAR spoofing attacks.

As shown in Table \ref{tb:svf_performance}, both SVF models achieve comparable AP compared to the original models. The marginal degradation comes from two-state training. More specifically, the distributions of the semantic segmentation outputs in the training and validation sets do not align well with each other. \CR{We find that the drop of AP will indeed cause a tiny amount of false negatives but will not influence the driving behaviors. Moreover, such degradation could be compensated with better training strategies (\textit{e.g.,}\xspace finer-tuning of the parameters) since the capacity of SVF is larger than the original models.}

\begin{table}[!h]
    \captionsetup{font=small}
    \footnotesize
    \caption{SVF-PointPillars' and SVF-PointRCNN's APs (\%) of 3D car detection on the KITTI validation set.}
    \vspace{-0.1in}
    \centering
    \begin{tabular}{|c|c|c|c|}
    \hline
    \multirow{2}*{Model} &  \multicolumn{3}{ c |}{Car Detection} \\
    \cline{2-4}
     &  Easy & Moderate & Hard \\
    \hline
    \hline
    SVF-PointPillars & 85.93 & 74.12 & 70.19 \\

    \hline
    SVF-PointRCNN &  88.12 & 76.56 & 74.81 \\

    \hline
    \end{tabular}
    \label{tb:svf_performance}
    \vspace{-0.3cm}
\end{table}

We then perform the robustness evaluation of SVF models. Figure \ref{fig:asr_svf} shows the $ASR$ of our proposed spoofing attacks. As shown, the attacks are no longer effective in SVF models. The $ASR$ reduces from more than 95\% (original models) to less than 4.5\% on both models with the maximum attack capability, which is also an around 2.2$\times$ improvement compared to CARLO-guarded models. The mean $ASR$ also drops from 80\% to around 2.3\%. We also perform ablation study on SVF, and demonstrate that the FV features are more important in SVF models (detailed in Appendix~\ref{ap:ablation}).

\begin{figure}[tp]
\vspace{-0.3cm}
\centering
\subfigure[$ASR$ of SVF-PointPillars.]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=0.95\linewidth]{images/usenix_svf_pointpillars.pdf}
\end{minipage}%
}%
\subfigure[$ASR$ of SVF-PointRCNN.]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=0.95\linewidth]{images/usenix_svf_pointrcnn.pdf}
\end{minipage}%
}%
\centering
\vspace{-0.4cm}
\caption{Attack success rates ($ASR$s) of proposed black-box spoofing attack on SVF models.}
\label{fig:asr_svf}
\vspace{-0.5cm}
\end{figure}

\nsubsubsection{Defense against White-box and Adaptive Attacks}
\label{ml_white_box}

Since SVF requires re-training for the model, we cannot directly evaluate Adv-LiDAR on SVF-Apollo (\S\ref{background_kitti}). As a result, we decouple the problem to whether Adv-LiDAR can fool both the semantic segmentation and 3D object detection modules. We first directly apply the attack traces that successfully fool Apollo 5.0 to the segmentation network and record the mean confidence score of all the points belonging to the attack trace. Figure \ref{fig:white_box_svf} shows that the mean confidence scores are consistently below 0.08, which is too low to be classified as a valid vehicle with mean confidence scores of around 0.73 in our trained model.

Model-level defenses are usually vulnerable to simple adaptive attacks~\cite{athalye2018obfuscated,carlini2017adversarial}. To demonstrate the effectiveness of SVF against adaptive attack, we assume that the adversaries are aware of the SVF architecture. The attack goal is to both fool the semantic segmentation and 3D object detection modules. We also leverage the formulation in \cite{cao2019adversarial} to utilize the global transformation matrix $H(\theta,\tau)$ to control the spoofed points.
\CR{
\begin{equation}
\begin{split}

    & \min_{\theta,\tau} \qquad - {\mathcal{L}}_{seg} (x \odot  V \cdot H(\theta,\tau)^{T}) \qquad

\end{split}
\label{eq:adaptive_2}
\end{equation}
where $\odot$ represents the point cloud merge in the front view and ${\mathcal{L}}_{seg}(\cdot)$ defines the average confidence score of the attack trace (\textit{i.e.,}\xspace $ V \cdot H(\theta,\tau)^{T}$).} Figure \ref{fig:adaptive_svf} shows that none of the attack traces' average confidence score reaches 0.2 in the segmentation module, which is still far from the mean average confidence score of valid vehicle 0.73. Therefore, the adaptive attacks also cannot break the robustness of SVF.

\section{Target Perception Models and Benchmark}
\label{Target}

\subsection{Target Perception Models}
\label{target_models}

\begin{figure}[thp!]
\centering
  \missingfigure[figwidth=8cm]{Testing a long text string}
  \caption{Comparison with State-of-the-art LiDAR-based Perception Architectures.}
  \label{fig:sota}
  \vspace{-0.5cm}
\end{figure}

\subsubsection{Baidu Apollo 5.0 \footnote{\url{https://github.com/ApolloAuto/apollo/tree/r5.0.0/modules/perception/lidar}}}

Figure.\ref{fig:sota} (a) shows the perception pipelines of Apollo model

\subsubsection{PointPillars \footnote{\url{https://github.com/nutonomy/second.pytorch}}} Figure.\ref{fig:sota} (b) shows the perception pipelines of PointPillars model

\subsubsection{PointRCNN \footnote{\url{https://github.com/sshaoshuai/PointRCNN}}}
Figure.\ref{fig:sota} (c) shows the perception pipelines of PointRCNN

\subsection{Target Benchmark}
\label{target_benchmark}

KITTI is

\subsection{Hypotheses Proposal}

\nsubsection{Experimental Validation}
\label{hypothesis_validation}

We design experiments (\textbf{E}) to test the existence of such potential erroneous predictions (\textit{i.e.,}\xspace \textbf{FP}) on three target models using the KITTI validation set.

\textbf{E1}: To validate \textbf{FP1}, we first randomly pick 100 point cloud samples $\mathcal{F}=\{F_{i}\}_{i=1}^{100}$ that contain 100 target \emph{occluded} vehicles $\{v_{i}\}_{i=1}^{100}$ with their point sets $V_i \subseteq F_i$. We then feed $\mathcal{F}$ into three target models and record the confidence scores (\textit{i.e.,}\xspace outputs of models to represent the confidence of detection) of the occluded vehicles as $s_i$ for each $v_i$.

Second, we leverage a global translation matrix $H(\theta,\tau)$ (Equation \ref{eq:trans}) to move every $V_{i}$ to a front-near location (\textit{i.e.,}\xspace 5-8 meters in front of the victim AV) in the point cloud $F_{i}$ as $V_{i}'$, where $\theta$ and $\tau$ correspond to the \emph{azimuth} and \emph{distance} of the translation, respectively:
\CR{
\vspace{-0.2cm}
\begin{equation}
\begin{split}
& \quad {V_{i}'}_{\wi} = {V_{i}}_{\wi} \\
\begin{bmatrix}
  {V_{i}'}_{\mathbf{w_x}}\\
  {V_{i}'}_{\mathbf{w_y}}\\
  {V_{i}'}_{\wz}\\
  1
\end{bmatrix} & = \begin{bmatrix}
  \cos{\theta} & -\sin{\theta} & 0 & \tau\cos{(\theta+\alpha)}\\
  \sin{\theta} & \cos{\theta} & 0 & \tau\sin{(\theta+\alpha)}\\
  0 & 0 & 1 & 0 \\
  0 & 0 & 0 & 1 \\
\end{bmatrix}
\cdot
\begin{bmatrix}
  {V_{i}}_{\mathbf{w_x}}\\
  {V_{i}}_{\mathbf{w_y}}\\
  {V_{i}}_{\wz}\\
  1
\end{bmatrix}
\end{split}
\label{eq:trans}
\end{equation}
}
$({V_{i}}_{\mathbf{w_x}},{V_{i}}_{\mathbf{w_y}},{V_{i}}_{\wz},{V_{i}}_{\wi})$ denotes the $xyz$-$i$ feature vectors (introduced in \S\ref{background_spoof}) of all points in $V_{i}$, and \CR{$\alpha=arctan ({V_{i}}_{\mathbf{w_y}} / {V_{i}}_{\mathbf{w_x}})$}. We make sure that there are no other objects standing between the LiDAR and each $V_{i}'$. By doing so, we construct a new set, where the points belong to the target occluded vehicles are moved to a front-near location by Equation \ref{eq:trans}. We further feed the new point cloud set $\mathcal{F}'$ into three target models and record the confidence scores of the translated points $V_{i}'$ as $s_{i}'$.

Experimental results show that all of the translated points $V_{i}'$ are detected by three target models, and we calculate the relative errors $e = \frac{|{s_i}'-{s_i}|}{{s_i}}$. Figure \ref{fig:cdf} shows the CDF of $e$ for three target models. As shown, 99.5\% of the picked occluded vehicles only have below 10\% fluctuations of their confidence scores, which successfully validate \textbf{FP1}.

The success of \textbf{E1} comes from the fact that LiDAR-based 3D object detection models perform \emph{amodal} perception, where given only the visible portions of a vehicle $v$, the model attempts to reason about occlusions and predict the bounding box for the complete vehicle (Figure \ref{fig:c1-3}). However, \CR{convolutional operations exploit spatial locality by enforcing a local connectivity pattern between neurons of adjacent layers. Such architecture thus ensures to produce the strongest response to a spatially local input pattern.} Since the \emph{occludee}'s and \emph{occluder}'s point sets $V$ and $\mathsf{O}(v)$ stand apart from each other in the 3D space, deep learning models may fail to identify the causality between $V$ and $\mathsf{O}(v)$ and thus learns to regress the bounding box for $v$ by $V$ only.

\textbf{E2}: To validate \textbf{FP2}, similarly, we first randomly pick 100 point cloud samples that contain 100 target \emph{distant} vehicles $\{v_{i}\}_{i=1}^{100}$ that locate farther than 30 meters away from the AV, and follow the same procedure with \textbf{E1} to record the confidence score changes. Experimental results show that all of the translated points $V_{i}'$ are detected by three target models, and we calculate the relative errors $e' = \frac{|{s_i}'-{s_i}|}{{s_i}}$. Figure \ref{fig:cdf} shows the CDF of $e'$ for three target models. As shown, 99.5\% of the picked distant vehicles only have below 7.5\% fluctuations of their confidence scores, which successfully validate \textbf{FP2}.

The success of \textbf{E2} comes from that 3D object detection models are designed to be non-sensitive to the locations of objects. For example, Apollo 5.0 does not incorporate location information in its hard-coded feature maps, and PointRCNN regards the centers of each bounding box as the origins of their coordinates. Hence the global locations of objects are not valued by the 3D object detection models in AD systems.

\begin{figure}[t!]
\centering
\vspace{-0.1cm}
  \includegraphics[width=1\linewidth, height=0.4\linewidth]{images/cdf.pdf}
  \vspace{-0.5cm}
  \caption{Left: CDF of the relative errors $e$ in \textbf{E1}. Right: CDF of the relative errors $e'$ in \textbf{E2}.}
  \label{fig:cdf}
  \vspace{-0.5cm}
\end{figure}

\nsubsection{Vulnerability Identification}
\label{vulnerability_id}

As mentioned earlier, the sensor attack capability $\mathcal{A}$ is far from spoofing a fully exposed front-near vehicle's point set. However, \textbf{E1} and \textbf{E2} provide two strategies for adversaries to launch spoofing attacks with fewer points and horizontal angles. As a result, attackers can directly spoof a vehicle imitating various occlusion and sparsity patterns that satisfy the sensor attack capability $\mathcal{A}$ to fool the state-of-the-art models. For example, the $V$ (red points) in Figure \ref{fig:c1-3} only contains 38 points and occupies 4.92$^\circ$ horizontally when translated to 6 meters in front of the AV. We confirm that it can deceive all three target models successfully, as visualized in Appendix~\ref{ap:figure}.

\begin{figure*}[!t]
\centering
\vspace{-0.25cm}
\begin{minipage}[b]{.72\linewidth}
\subfigure[$ASR$ of Apollo 5.0.]{
\centering
\includegraphics[width=0.3333\linewidth]{images/usenix_apollo_asr.pdf}
}%
\subfigure[$ASR$ of PointPillars.]{
\centering
\includegraphics[width=0.3333\linewidth]{images/usenix_pointpillars_asr.pdf}
}%
\subfigure[$ASR$ of PointRCNN.]{
\centering
\includegraphics[width=0.3333\linewidth]{images/usenix_pointrcnn_asr.pdf}
}%
\end{minipage}
\hfill
\begin{minipage}[b]{.25\linewidth}
	\centering
	\includegraphics[width=\linewidth]{images/mesh.pdf}

\end{minipage}

\centering
\vspace{-0.2cm}

\begin{minipage}[b]{.70\linewidth}
\caption{Attack success rates ($ASR$s) of proposed black-box spoofing attacks on target state-of-the-art models.}
\label{fig:asr_all}
\end{minipage}
\hfill
\begin{minipage}[b]{.26\linewidth}
\centering
\caption{The process of generating attack traces for $\mathcal{R}$ from the implemented renderer.}
\label{fig:render}
\end{minipage}
\vspace{-0.8cm}
\end{figure*}

The vulnerability comes from the observation that the state-of-the-art 3D object detection architectures ignore the distinct physical features of LiDAR. Therefore, they leave a gap, as well as an attack surface, between the model capacity and LiDAR point clouds. We further abstract the neglected physical features as two occlusion patterns inside the LiDAR point clouds, described below.

\textbf{Inter-occlusion}. We abstract the typical occlusion introduced in \S\ref{hypothesis} as inter-occlusion. As its name indicates, inter-occlusion describes a causal relationship between \emph{occludee} and the corresponding \emph{occluders} (\textit{i.e.,}\xspace the \emph{occluders} cause the \emph{occludee} partially visible). \textbf{FP1} violates the physical law of inter-occlusion since a translated ``occluded'' vehicle's point set $V'$ no longer has its valid \emph{occluder} $\mathsf{O}(v)$. However, \textbf{E1} demonstrates that state-of-the-art LiDAR-based perception models overlook such inter-occlusions in the point clouds.

\textbf{Intra-occlusion}. We abstract the other occlusion pattern hidden inside an object as intra-occlusion. The facing surface of a solid object (\textit{e.g.,}\xspace~a vehicle) occludes itself in the point cloud, which indicates that the LiDAR cannot perceive the interior of the object (Figure \ref{fig:fsd}). \textbf{FP2} violates the physical law of intra-occlusion since the abnormal sparseness of a translated ``distant'' vehicle's point set $V'$ can no longer fully occlude a valid vehicle since other laser pulses could penetrate its ``surface''. However, \textbf{E2} demonstrates that state-of-the-art LiDAR-based perception models \CR{are unable to differentiate reflected points of real solid objects from sparse injected points of the same overall shape so that they also} overlook the intra-occlusions in the LiDAR point clouds.

To demonstrate the potential real-world impacts of this identified vulnerability, we construct the first black-box spoofing attack on state-of-the-art LiDAR-based perception models in \S\ref{bb_attack}. We find the violations of the physical law of occlusion generally enable LiDAR spoofing attacks. Therefore, we perform the first defense study, exploiting the occlusion patterns as physical invariants to detect spoofing attacks in \S\ref{piad}. Lastly, in \S\ref{ml_defense}, we present a general architecture for robust LiDAR-based perception that embeds occlusion patterns as robust features into end-to-end learning.

\section{A General Design-level Vulnerability}
\label{vulnerability}

As introduced in \S\ref{background_kitti}, a point cloud sample in KITTI consists of around 133,000 points, where a front-near vehicle should contain at least 2000 points and occupy about 15 $^\circ$ of horizontal angles\cite{Geiger2013IJRR}. However, Cao \textit{et al.}\xspace demonstrated that an attack trace with 60 points and 8 $^\circ$ of horizontal angles is sufficient to spoof a front-near vehicle in Apollo 2.5\cite{cao2019adversarial}. It still remains unclear why such spoofing attacks can succeed even with a huge gap existing between adversarial and normal vehicles. Despite the promising results of attacking image-based machine learning models shown in previous studies\cite{eykholt2018robust,220580}, their approaches have internal differences with Adv-LiDAR. We distinguish the spoofing attacks towards LiDARs with extensive prior work in mainly three aspects:

\textbf{Different perturbation methods.} Images and point clouds have different data structures, which further leads to different perturbation methods applied. Images have compact structures (\textit{i.e.,}\xspace $H \times W \times C$), where $H \times W$ stands for the image size and $C$ contains the RGB information. However, point clouds are unordered which are oftentimes represented as $N \times C$, where $N$ is the number of points and $C$ contains the location and intensity information (\textit{i.e.,}\xspace$xyz$-$i$). Attackers are able to generate adversarial examples by modifying the RGB values in images. For LiDAR-based attacks, however, attackers can directly shift the point in the 3D Euclidean space as long as it obeys the physics of LiDAR (\S\ref{background_direct_spoof}).

\textbf{Different perturbation capabilities.} Prior attacks towards 2D images always treat the whole target area as the attack surface since the threat models assumes that attackers have full controls over the target object (\textit{e.g.,}\xspace the attackers can potentially modify any areas of a stop sign\cite{eykholt2018robust}). However, due to the characteristics of LiDAR spoofing attacks, the attack surface is limited by the sensor attack capability ($\mathcal{A}$) in \S\ref{threat}. Such a small attack surface also introduces more difficulties to launching the adversarial spoofing attack.

\textbf{Different perturbation constraints.} Prior attacks towards 2D images leverage \emph{l-p} norms as main constraints for the formulated optimization problem\cite{eykholt2018robust} whose goal is to minimize the perturbation to be stealthy. Such constraints are not applicable to attacking LiDAR because point clouds are not perceived by human beings, thus, stealthiness is not a focus in such attacks. However, due to the limited sensor attack capability (\S\ref{threat}), the optimized attack traces must not exceed the capability boundary, in which case the sensor attack capability naturally becomes the main constraint for attacking LiDAR.

Therefore, in order to unveil the general vulnerability exposed by Adv-LiDAR, it is necessary to take the distinct characteristics of the LiDAR point cloud into considerations, which are different from image-based attacks.

\subsection{Hypotheses Proposal}
\label{hypothesis}

\begin{figure}[thp!]
\centering
  \includegraphics[width=6cm]{images/C1-3.pdf}
  \caption{Illustration of \textbf{C1-3} of LiDAR point cloud.}
  \label{fig:c1-3}
\end{figure}

Although deep learning has been successfully applied to LiDAR-based object detection in autonomous driving, the interpretability of such models is always a bottleneck factor. In this section, we first provide several empirical observations (\textbf{O}) of LiDAR-based perception models and unique characteristics (\textbf{C}) of LiDAR point cloud data, based on what, we further propose two hypotheses (\textbf{H}) that could lead to a design-level vulnerability of state-of-the-art models.

\textbf{C1}: Occlusion between objects is inevitable in LiDAR point clouds and will cause partial visibility of the occluded objects. As introduced in \S\ref{background_spoof}, a LiDAR sensor functions by firing laser pulses and capturing their returns. As a result, each point in a point cloud represents the distance to the \emph{nearest} solid object that along the laser ray. Similar to human eyes, a LiDAR can only perceive parts of an object (\textit{e.g.,}\xspace a vehicle) if other obstacles standing between the LiDAR and the object that exclude the laser beams.

In this paper, we describe such occluded objects as \emph{occludees} and the objects that occlude others as \emph{occluders}. Particularly, we use $\mathsf{O}(x)$ to express the set of points that occludes the object $x$, and we use the captain letter $X$ to express the set of points that belong to object $x$ in the point cloud $F$ ($X = \{\vec{p}\, | \,\vec{p} \in x\} \subseteq F$, where $\vec{p} \in \mathbb{R}^3$).

\alfred{need a figure to illustrate the point} \mytodo{I will draw a figure to show C1-3}

\textbf{O1}: Deep learning identifies correlations in the training data, but it arguably cannot build causal models by purely observing correlations\cite{doshi2017towards}. The state-of-the-art LiDAR-based perception models (\S\ref{background}) most likely preserve this characteristic as they consist of basic fully-connected and convolutional layers.

\textbf{C2}: The \emph{occludee} and the corresponding \emph{occluders} described in \textbf{C1} stand separately from each other in the 3D Euclidean space. Intuitively, an occlusion relation exists as long as the \emph{occluders} shades the \emph{occludee} from the perspective of the LiDAR sensor As a result, $X$ and $\mathsf{O}(x)$ may be far apart with each other in the point cloud.

\textbf{O2}: CNNs exploit spatial locality by enforcing a local connectivity pattern between neurons of adjacent layers. The architecture thus ensures to produce the strongest response to a spatially local input pattern. PointPillars and Apollo 5.0 both leverage CNNs for the final detection task, so they are supposed to remain such characteristic, too. PointRCNN only accepts points inside region proposals into its second stage pipeline, thus spontaneously preserves locality.

Occlusion is a causation relationship. The existence of the \emph{occluders} leads to the partial visibility of the corresponding \emph{occludee}. From \textbf{O1}, deep learning models can barely identify such causation relationship introduced by occlusions. Besides, since the $X$ and $\mathsf{O}(x)$ stand independently from each other in the 3D Euclidean space (\textbf{C2}), LiDAR-based perception models have limitations to even identify the correlation between the $X$ and $\mathsf{O}(x)$ due to \textbf{O2}. Thus, we propose our first hypothesis:

\vspace{0.2cm}

\textbf{H1}: \emph{Given a point cloud $F$ containing an occluded object $x$, if $\mathcal{M}(F) \Rightarrow \{\mathsf{B}(x),\ldots\}$, we hypothesize $\mathcal{M}(F-\mathsf{O}(x)) \Rightarrow \{\mathsf{B}(x),\ldots\}$, where $\mathsf{B}(x)$ denotes the bounding box parameters of $x$.}

\vspace{0.2cm}

\textbf{C3}: Due to the working principles of LiDAR sensors (\S\ref{spoof}), the generated point clouds are collected uniformly in vertical and horizontal angles but not in Euclidean space. The points farther away from the LiDAR sensor are generally sparser than the points nearby the sensor. The horizontal and vertical resolutions are  $\propto r\times\sin{\Delta \theta}$ and $\propto r\times\sin{\Delta \phi}$ respectively, where $r$ is the distance from the point to the LiDAR sensor, and $\Delta \theta$ and $\Delta \phi$ are the horizontal and vertical fire angle intervals.

\textbf{O3:} CNN-based 2D object detection models are oftentimes translation equivariant\cite{ren2015faster}. Intuitively, if an object is translated to another location in the image plane, its detected bounding box should shift with the translation of the original object\cite{weiler20183d}:
\begin{equation}
    \mathcal{M} \circ \mathsf{T}\,_{\vec{t}} = \mathsf{T}\,_{\vec{t}} \circ \mathcal{M}
\end{equation}
where $\mathsf{T}\,_{\vec{t}}(\vec{x}) := (\vec{x} + \vec{t})$, where $\vec{x},\vec{t} \in \mathbb{R}^2$ and $\mathcal{M}$ is the deep learning model. In terms of LiDAR-based perception models, first, the hard-coded features in the BEV-based models (\textit{e.g.,}\xspace Apollo 5.0) do not contain location information in $xy$ plane, therefore, will most likely possess the characteristic in $xy$ plane. Second, for the voxel-based models, the initial features inside each voxel contain the location information (\textit{i.e.,}\xspace $x,y,z$). However, they also incorporate the local location $x_c,y_c,z_c$ where the $c$ subscript denotes distance to the arithmetic mean of all points in the voxel (\S\ref{background}) which are translation equivariant. Thus, the characteristic may still hold in the voxel-based models. Third, in the second stage of the point-wise models, canonical coordination transformations are performed in each region proposal, regarding the centers of the region proposals as the origins of their coordinates. This canonical transformation automatically makes the model translation equivariant.

Even though the density of point clouds varies a lot along their distances to the sensor (\textbf{C3}), we propose our second hypothesis based on \textbf{O3}:

\vspace{0.2cm}

\textbf{H2}: \emph{Given a point cloud $F$ containing an object $x$, if $\mathcal{M}(F) \Rightarrow \{\mathsf{B}(x),\ldots\}$, we hypothesize $\mathcal{M}((F - X) \cup \mathsf{T}\,_{\vec{t}}(X)) \Rightarrow \{\mathsf{B}(\mathsf{T}\,_{\vec{t}}(x)),\ldots\}$, where $\mathsf{B}(x)$ denotes the bounding box parameters of $x$.}

\subsection{Hypotheses Validation}
\label{hypothesis_validation}

We test the two proposed hypotheses using point clouds from the KITTI validation set.

\textbf{V1}: To validate \textbf{H1}, we first randomly pick 100 point cloud frames $\mathcal{F}=\{F_{i}\}_{i=1}^{100}$ that contain 100 target \emph{occluded} vehicles $\{v_{i}\}_{i=1}^{100}$ where the corresponding set of points $V_i \subseteq F_i$ respectively. We then feed $\mathcal{F}$ into three target models and record the confidence scores of the occluded vehicles as $s_i$ for each $v_i$.

Second, we eliminate the specific $\mathsf{O}(v_i)$ for each vehicle $v_i$. The rest of points of each point cloud $F_{i}$ form a new set of point cloud frames $\mathcal{F}^{'}=\{F_{i}-\mathsf{O}(v_{i})\}_{i=1}^{100}$. We feed $\mathcal{F}^{'}$ into the three target models and record new confidence scores of the occluded vehicles as ${s_i}'$ for each $v_i$.

Figure \ref{fig:cdf} shows the CDF of $e = \frac{|{s_i}'-{s_i}|}{{s_i}}$ for all three target models. As shown, 99.5\% of the picked occluded vehicles only have below 2\% fluctuations of their confidence scores when their corresponding $\mathsf{O}(\cdot)$ are eliminated from the point cloud on all target models. The experimental results validate \textbf{H1}.

\textbf{V2}: To validate \textbf{H2}, we first generate a point cloud $G$ simulating the ground plane's points as reference for the $xy$ plane. Then we again randomly pick and extract 100 vehicles' point clouds $\mathcal{V}=\{V_{i}\}_{i=1}^{100}$. It is worth noting that they are not necessary occluded. We further merge the point cloud of the ground plane $G$ with $\mathcal{V}$ forming a new set $\mathcal{F}=\{G \cup V_{i}\}_{i=1}^{100}$. We feed $\mathcal{F}$ into three target models and record the confidence scores of the occluded vehicles as $s_i$ for each $v_i$.

Second, we grid the ground plane $G$ into 1 $m^2$ cells $\{(x_{p},y_{q})\}_{p,q=1}^{N_{x},N_{y}}$, where $N_{x}$ and $N_{y}$ are the number of cells in $x$ and $y$ axes respectively. By applying a global transformation matrix $H(\theta,\tau)$\cite{cao2019adversarial} to $\mathcal{V}$, we translate every $V_{i}$ to every cell $(x_{p},y_{q})$ on the ground plane as $V_{i,p,q}'$:
\begin{equation}
\begin{split}
& \quad {V_{i,p,q}'}_{\wi} = {V_{i}}_{\wi} \\
\begin{bmatrix}
  {V_{i,p,q}'}_{\mathbf{w_x}}\\
  {V_{i,p,q}'}_{\mathbf{w_y}}\\
  {V_{i,p,q}'}_{\wz}\\
  1
\end{bmatrix} & = \begin{bmatrix}
  \cos{\theta} & -\sin{\theta} & 0 & \tau\cos{\theta}\\
  \sin{\theta} & \cos{\theta} & 0 & \tau\sin{\theta}\\
  0 & 0 & 1 & 0 \\
  0 & 0 & 0 & 1 \\
\end{bmatrix}
\cdot
\begin{bmatrix}
  {V_{i}}_{\mathbf{w_x}}\\
  {V_{i}}_{\mathbf{w_y}}\\
  {V_{i}}_{\wz}\\
  1
\end{bmatrix}
\end{split}
\end{equation}
where $({V_{i}}_{\mathbf{w_x}},{V_{i}}_{\mathbf{w_y}},{V_{i}}_{\wz},{V_{i}}_{\wi})$ denotes the $xyz$-$i$ features of all points in $V_{i}$. We merge the $V_{i,p,q}'$ with $G$ forming $\mathcal{F}'=\{V_{i,p,q}' \cup G\}_{i,p,q=1}^{100,N_{x},N_{y}}$. We feed $\mathcal{F}'$ into three target models and record the confidence score of the translated vehicle as $s_{i,p,q}'$.

Figure \ref{fig:cdf} also shows the CDF of $e' = \frac{|s_{i,p,q}'-{s_i}|}{{s_i}}$ for all three target models. As shown, the relative error $e'$ are larger than $e$, which is expected since the models contain location-sensitive features (\textbf{O3}). However, the results the fluctuations are still below 10\% for 99.5\% picked vehicles, and we verify that with such fluctuations, the models still detect these vehicles as valid which validate \textbf{H2}.

\begin{figure}[thp!]
\centering
  \includegraphics[width=1\linewidth, height=0.4\linewidth]{images/cdf.pdf}
  \caption{CDF of the relative errors $e$ (left) and $e'$ (right).}
  \label{fig:cdf}
  \vspace{-0.5cm}
\end{figure}

\subsection{Vulnerability Identification}
\label{vulnerability_id}

In this section, we unveil how the validated hypotheses (\textbf{V1-2}) lead to a design-level vulnerability of state-of-the-art LiDAR-based perception models.

As introduced in \S\ref{background_adv_spoof}, Equation \ref{eq:adv_lidar} formulates the optimization problem of the adversarial spoofing attack, where $\mathcal{A}$ is the sensor attack capability.

Instead of spoofing a complete vehicle, to spoof an occluded vehicle is a much easier task since an occluded vehicle occupies smaller horizontal angles and less points than a complete one. \textbf{V1} shows that even without a valid \emph{occluder}, an ``occluded'' vehicle can still be detected by the models. Moreover, attackers can spoof even sparser points of an occluded vehicle to further reduce the cost. Due to the non-uniform density of LiDAR point clouds (\textbf{C3}), a vehicle far away from the LiDAR sensor contains sparse points. \textbf{V2} shows that translating such a distant vehicle's points to a nearby location will make the models detect the translated points as a vehicle at the new location.

As a result, instead of formulating the $\oplus$ function and solving Equation \ref{eq:adv_lidar} by gradient descent-based methods, attackers can directly spoof an ``occluded'' vehicle with sparse points that satisfy the sensor attack capability $\mathcal{A}$ to fool the state-of-the-art models. \mytodo{not sure whether the explanation is clear.}