# Time-reversed photoacoustic guided time-reversed ultrasonically encoded
optical focusing

Juze Zhang  Hybrid Imaging System Laboratory, Shanghai Engineering Research
Center of Intelligent Vision and Imaging, School of Information Science and
Technology, ShanghaiTech University, Shanghai 201210, China  Zijian Gao
Hybrid Imaging System Laboratory, Shanghai Engineering Research Center of
Intelligent Vision and Imaging, School of Information Science and Technology,
ShanghaiTech University, Shanghai 201210, China  Xiaopeng Yu  Hybrid Imaging
System Laboratory, Shanghai Engineering Research Center of Intelligent Vision
and Imaging, School of Information Science and Technology, ShanghaiTech
University, Shanghai 201210, China  Peng Ge  Hybrid Imaging System Laboratory,
Shanghai Engineering Research Center of Intelligent Vision and Imaging, School
of Information Science and Technology, ShanghaiTech University, Shanghai
201210, China  Feng Gao  Hybrid Imaging System Laboratory, Shanghai
Engineering Research Center of Intelligent Vision and Imaging, School of
Information Science and Technology, ShanghaiTech University, Shanghai 201210,
China  Fei Gao  Hybrid Imaging System Laboratory, Shanghai Engineering
Research Center of Intelligent Vision and Imaging, School of Information
Science and Technology, ShanghaiTech University, Shanghai 201210, China
Corresponding author<EMAIL_ADDRESS>
###### Abstract

Deep-tissue optical imaging is a longstanding challenge limited by scattering.
Both optical imaging and treatment can benefit from focusing light in deep
tissue beyond one transport mean free path. Wavefront shaping based on time-
reversed ultrasonically encoded (TRUE) optical focusing utilizes ultrasound
focus, which is much less scattered than light in biological tissues as the
’guide star’. However, the traditional TRUE is limited by the ultrasound
focusing area and pressure tagging efficiency, especially in acoustically
heterogeneous medium. Even the improved version of iterative TRUE comes at a
large time consumption, which limits the application of TRUE. To address this
problem, we proposed a method called time-reversed photoacoustic wave guided
time-reversed ultrasonically encoded (TRPA-TRUE) optical focusing by
integrating ultrasonic focus guided by time-reversing PA signals, and the
ultrasound modulation of diffused coherent light with optical phase
conjugation (OPC), achieving dynamic focusing of light into scattering medium.
Simulation results show that the focusing accuracy of the proposed method has
been significantly improved compared with conventional TRUE, which is more
suitable for practical applications that suffers severe acoustic distortion,
e.g. transcranial optical focusing.

††journal: optica

##  1 Introduction

Optical imaging in deep tissues is an open challenge in biomedical imaging
research. The focus of light in deep tissues can not only enhance the output
signal intensity of optical imaging systems (such as fluorescence imaging),
but also improve imaging sensitivity and depth. On the other hand, it can also
be directly applied to optogenetics, photodynamic therapy, laser ablation and
other fields that require light energy concentration [1, 2, 3].

However, one longstanding challenge in optical imaging is scattering. Strong
scattering causes the wavefront to be distorted as they pass through the
tissue. Meanwhile, the photons are absorbed by the tissue during propagation,
which greatly limits the penetration of light into deeper biological tissues.
When a beam of light enters the tissue, photons will bounce many times, and
its movement will become random after transport mean free path [4]. Until
recently, most optical techniques used for imaging inside tissues mainly rely
on only un-scattered (i.e. ballistic) photons such as confocal microscopy,
two-photon microscopy and optical coherence tomography (OCT) [5]. Therefore,
the ability to focus light inside the scattering medium can revolutionize the
deep-tissue non-invasive fluorescence microscope, optical tweezers,
optogenetics, microsurgery, and phototherapy, covering many areas of
biophotonic applications[6].

Wavefront shaping may be the best method to address this problem by shaping
the wavefront of the incident light to counteract the scattering [7]. Three
main classes of wavefront shaping approaches are being actively developed,
including model-based wavefront shaping [8], feedback-based wavefront shaping
[9] and optical time reversal/optical phase conjugation (OPC) [10, 11].

Model-based wavefront shaping depends on the numerical computation using a
digital model of the sample. The microscopic refractive index model of the
sample, which serves as the input for the calculations, is obtained from the
image data itself [8]. Feedback-based wavefront shaping depends on the
detection of the intensity feedback [12] at a desired location to find the
optimal wavefront maximizing that feedback signal, by using some advanced
optimization algorithm such as genetic algorithm [12], and even deep learning
[13]. This feedback signal is often obtained by using a camera to measure the
scattered light, but it is unrealistic to place a camera inside biological
tissue. Lai et al [14]. proposed to use the photoacoustic signal as a feedback
signal to optimize the wavefront of light, which makes the feedback-based
method suitable for biological tissues. Despite this, the process of focusing
consumes a lot of time to optimize, which limits its wide applications.

Figure 1:  Simplified schematic diagram of principle for TRUE optical
focusing. (a) Initially, a wide laser beam or a narrow laser beam is incident
on the scattering medium. After entering the scattering medium, billions of
photons undergo some scattering events, and the path of the photons becomes
random. When the scattered photon passes through the position under
insonification, part of the photons are tagged by the ultrasound, resulting in
a frequency shift. The white solid line represents the photon path of the
tagged part, and finally reaches the SLM and then is recorded. (b) Through the
holographic method, we can calculate the phase of the photon in the tagged
part, take the conjugate of the phase and load it on the SLM. As a result, the
photons will follow their path and return to the region under insonification,
just like going back in time.

Alternatively, in OPC, the optimal wavefront is obtained from a single
measurement of the scattered photons propagating from a source located behind
or inside the turbid medium. Subsequently, a focus is formed by playing back
the conjugate of this field using a phase conjugate mirror. Similar to
feedback-based wavefront shaping, OPC requests a ’guide star’ in the
scattering medium to guide the wavefront focusing [15, 16]. To address this,
Xu et al [17] proposed a scheme termed time-reversed ultrasonically encoded
(TRUE) optical focusing, which uses ultrasound focus as the ’guide star’.
Figure 1 illustrates the process of TRUE. Unlike light, ultrasound is much
less scattered in biological tissues, which can create a virtual source of
light frequency shifted by the acousto-optic effect as the ’guide star’. The
ultrasound emitted by a focused ultrasound transducer modulates scattered
photons (ultrasonically encoded) (Fig. 1a). As a result, the encoded photons
will have a frequency shift to form a ’guide star’ in the ultrasound focus
area. As a ’new’ light source, the ’guide star’ forms an interference pattern
with the reference light and is recorded by the holographic recording medium.
When the hologram is reproduced, the phase-conjugated object light will return
to the ultrasonic modulation area through time reversal, forming an optical
focus inside the scattering medium (Fig. 1b).

The resolution of the early TRUE technology is determined by the size of the
ultrasound focus (25-50$\mu$m) [15]. However, the focusing ability of
ultrasound is essential to ensure sufficient localized energy deposition and
sufficient number of photons modulated in the target area (see Section 2,
Methods). When the target area is located behind an acoustically heterogeneous
layer (such as an irregular layer of adipose tissue, or bone), the acoustic
focusing performance becomes the main bottleneck limiting its wide
application. One possible solution that has been proposed to solve the TRUE
resolution problem is iterative TRUE [18, 19, 20], that is, the feedback loop
between the conjugate field and the detection field can reduce the focus size
by three times validated through experiments. Nevertheless, the time consuming
of iterative schemes will increase the complexity of the system.

Recently, photoacoustic (PA) generation in the target have been proposed as an
alternative, e.g. Fink et al. have successfully verified the feasibility of
time-reversal of PA signals [21]. Moreover, our previous works have also
demonstrated the potential of ultrasonically focus guided by time-reversing
the transcranial PA signals of an optically selective target in a nonselective
background [22].

In this paper, we propose a method by hybrid integrating photoacoustic,
ultrasonic time-reversal, and OPC to achieve real-time ultrafast 2D/3D optical
focusing in both optically and acoustically heterogeneous medium with one
snapshot.

##  2 Method

Our proposed method, called time-reversed photoacoustic wave guide time-
reversed ultrasonically encoded (TRPA-TRUE) optical focusing, integrates
ultrasonic focus guided by time-reversing PA signals, and the ultrasonic
modulation of diffused coherent light with OPC, achieving dynamic focusing of
light into a scattering medium. Figure 2 illustrates the concept of our
method. First, an unfocused light from a laser source ( $\lambda$=532 nm) is
illuminated into the medium. Optically absorbing regions act as acoustic
sources through the photoacoustic effect, which generates a photoacoustic wave
detected by transducer array (see Fig. 2a). Subsequently, we perform acoustic
time reversal operation on the received PA signals to achieve ultrasound
focusing inside the medium (Fig. 2b). When the time-reversed PA signal reaches
the optical absorption regions, part of the photons will be tagged in this
snapshot in the optical absorption area and then reach the OPC (Fig. 2c).
Finally, we perform TURE operation to complete whole process of TRPA-TRUE,
which achieves real-time ultrafast 2D/3D optical focusing in both optically
and acoustically heterogeneous medium with one snapshot.

In the following sections, we will start with the acousto-optic effect and
explain why the focus of ultrasound determines our TRUE effect, then we will
explain the feasibility of time reversal of the PA signal. Finally, we will
also perform a simulation study to verify our method.

Figure 2:  Simplified schematic diagram of principle for TRPA-TRUE optical
focusing. (a) First, an unfocused light from a laser source ( $\lambda$ =532
nm) is illuminated into the medium and optically absorbing regions act as
acoustic sources through the photoacoustic effect, which generates
photoacoustic wave and then be detected. (b) Subsequently, we perform time
reversal operation on the received PA signal to achieve ultrasound focusing.
(c) When the time-reversed PA signal reaches the optical absorption regions, a
part of the photons will be tagged in this snapshot in the optical absorption
area and then reach the OPC. (d) Finally, we perform TURE operation to
complete TRPA-TURE.

###  2.1 Acousto-Optic Interaction within Turbid Media

There are currently three widely recognized modes to interpret the AO effect:
The displacement of the scatterers in the turbid medium caused by the acoustic
field [23], the disturbance of the refractive index of the medium caused by
the mechanical strain due to the acoustic field [24], and perturbation of the
optical properties of the medium due to compression and rarefaction of the
medium by the acoustic field [25]. However, the last mechanism is very weak,
so we can ignore it here [26].

####  2.1.1 Ultrasound-induced change of the refractive index

Denote the acoustic pressure field at position $\bf{r}$ and time $t$ as the
following equation:

$P(\mathbf{r},t)=P_{0}\sin(\omega_{a}t-\mathbf{k_{a}}\cdot\mathbf{r})$ (1)

where $P_{0}$ is the acoustic amplitude, $\omega_{a}$ is the acoustic angular
frequency and $\mathbf{k_{a}}$ is the acoustic wave vector.

Acoustic waves propagating through the medium produce strain, which alters the
local permittivity of medium. The adiabatic piezo-optic coefficient
$\eta=\frac{\partial n}{\partial p}$ relates the change in the index of
refraction $n$ to the acoustic pressure $P(\mathbf{r},t)$. The refractive
index change can be described as:

$\Delta n=n_{0}(\frac{\partial n}{\partial p})P(\mathbf{r},t)$ (2)

where $n_{0}$ is the background refractive index without ultrasound
modulation.

####  2.1.2 Ultrasound-induced displacement of scatterers.

To describe the locations of any particles, we model the ultrasound-induced
displacement of particle at position $\bf{r}$ and time $t$ by using the
following equation:

$A(\mathbf{r},t)=A_{0}\mathbf{\hat{k}_{a}}\sin(\omega_{a}t-\mathbf{k_{a}}\cdot\mathbf{r}+\Phi)$
(3)

where $A_{0}$ is the displacement amplitude, $\mathbf{\hat{k}_{a}}$ is the
unit vector of $\mathbf{k_{a}}$, and $\Phi$ is particular phase offset between
particle displacement and ultrasound pressure. Recently, Huang et al. have
investigated the ultrasound–light interaction in scattering media[27], and we
can know from their theory that $A_{0}=\frac{P_{0}}{\omega_{a}\rho v_{a}}$,
where $\rho$ is the density of the background medium and $v_{a}$ is the
velocity of acoustic wave. In their theory, phase offset $\Phi$ can be
approximately replaced by $3\pi/2$ between particle displacement and
ultrasound pressure when when particles move.

Combining (2) and (3), we can easily get that when the Monte Carlo(MC)
simulation is performed, the number of the ultrasound encoded photons shows
positive correlation with acoustic intensity, whether due to ultrasound-
induced change of the refractive index or ultrasound-induced displacement.

####  2.1.3 Ultrasound tagging efficiency in scattering media

Another important factor, tagging efficiency, should be taken into account.
Here we use $\eta$ to represent tagging efficiency and we give its
mathematical expression as following:

$\eta=\frac{P_{\text{tag}}}{P_{\text{untag}}+P_{\text{tag}}}$ (4)

where $P_{\text{tag}}$ is the power of light that is tagged by ultrasound
pressure and $P_{\text{untag}}+P_{\text{tag}}$ is the total power of light
that pass through the ultrasound pressure field. It is worth noting that here
we only consider the part of the light passing through the ultrasound pressure
field, not all the photon power. Unfortunately, we did not derive the
mathematical expression between the tagging efficiency and the ultrasound
pressure field, so here we prove the relationship between them through
simulation.

To calculate the tagging efficiency, we only need to consider power spectrum
of one speckle exiting the skull and we represent it by electric field $E(t)$.
We have $E(t)=\sum_{i}E(\theta_{i},t)$, where $\theta_{i}$ represents the
$i$’th photon in this speckle. Because in MC simulation, each photon is
independent, we did not treat the photon as a wave phenomenon and ignore
features such as phase and polarization. So we can calculate power spectrum by
the following equation[27]:

$\mathbb{E}\\{|\mathcal{F}[E(t)]|^{2}\\}=\mathbb{E}\\{\sum_{i}|\mathcal{F}[E(\theta_{i},t)]|^{2}\\}$
(5)

where $\mathbb{E}$ represents expectation and we can calculate power spectrum
by $|\mathcal{F}[E(t)]|^{2}$. With the above definition, we can estimate the
tagging efficiency from the overall average of the power spectrum of all
speckles. But calculating all speckles seems to be too computationally
expensive, so we can estimate the tagging efficiency by sampling a subset of
photon through the scattering medium.

Thanks to Yujia Huang for opening up their ultrasound-tagging code on github
and uploading their experimental data[27]. We downloaded the code and tested
it to get the relationship between the tagging efficiency and ultrasound
pressure as well as the ultrasound frequency. The comparison result is shown
in Fig. 3.

Figure 3:  Tagging efficiency with 1 MHz, 2.25MHz and 3.5MHz acoustic
frequencies versus acoustic pressure. Figure 4:  The diagram of the simulation
environment modeling process. (a) Get each slice from MRI data. (b) The MRI
data of brain were imported into MATLAB and segmented into air, water,
midbrain, white matter, grey matter, spinal fluid scalp and skull using
thresholds segmentation. (c) After segmenting each slice, we stack the
processed MRI images to form a 3-D model, and then load it into the k-wave
simulation package.

Figure 3 shows that the higher the ultrasound pressure, the higher the
ultrasound tagging efficiency. This is also consistent with our prediction
that the ultrasound focus determines the TRUE focusing ability. And at the
same time, we can also get that as the frequency increases, the tagging
efficiency is decreasing. According to what we have learned, in order to
obtain a smaller ultrasound focus area, many research groups use a higher
ultrasound frequency (e.g. 50 MHz), but this will result in a decrease in the
tagging efficiency of ultrasound, which increases the difficulty of detection.
In addition, even if we ignore the negative effects of the reduction of
ultrasound tagging, the use of higher ultrasound frequencies greatly reduces
penetration depth, because high-frequency acoustic signals tend to be more
easily attenuated by tissues. To sum up the above points, in order to achieve
better optical focusing through the TRUE system, we first need to improve the
efficiency of the ultrasound tagging ‘guide star’, and our method can greatly
improve this performance.

###  2.2 Time invariance of PA signals

Optically absorbing regions act as acoustic sources through the photoacoustic
effect described by the following equation[28]

$\Delta p-\frac{1}{c^{2}}\frac{\partial^{2}p}{\partial
t^{2}}=-\frac{\beta}{c_{p}}\frac{\partial H}{\partial t}$ (6)

where $p$ is the acoustic pressure, $c$ is the speed of sound, $\beta$ is the
isobaric thermal expansion coefficient, $c_{p}$ is the specific heat, and $H$
is the volumetric density of optical energy per unit time. If $p(r,t)$ is a
solution of the equation, then $p(r,-t)$ is also a solution. In other words,
the equation is invariant under time reversal[29]. Thus, diverging waves
emanating from a source can be reversed in time and back-propagated in space
to focus at the source.

##  3 Simulation of TRPA-TRUE on a Brain model

Our method was verified by performing simulations of TRPA-TRUE within Turbid
Media using a brain model (Fig. 4a), which is downloaded from the Scalable
Brain Atlas [30]. The acoustic simulations are implemented in k-space using
the k-Wave toolbox of MATLAB (R2019a, Mathworks Inc., Natick, MA, USA) [31],
which provides raw PA data, and simulates the process of time-reversed sound
wave propagation. In the optical simulation part, we use the Monte Carlo (MC)
photon propagation model [32], which has been widely used in acousto-optic
imaging (AOI) simulation [27, 33, 34], because it can directly simulate the
interaction of acousto-optic in the medium. The simulation medium proprieties
can be set using masked operations assigning the values from the brain model.
The properties of the different intracranial propagation medium is shown in
Table 1.

As shown in Fig. 4, the MRI data of brain were imported into MATLAB and
segmented into air, water, midbrain, white matter, grey matter, spinal fluid
scalp and skull using thresholds segmentation. After segmenting each slice, we
stack the processed CT images to form a 3-D model, and then load it into the
k-wave simulation package. Figure 4 illustrates the steps involved in the
segmentation and acoustic brain simulation configuration process. The
simulation setup in Fig. 4(c) shows that 64 transducers were placed around the
center of $xz$ plane to acquire the PA signals from the sample. The distance
from the center to the acoustic sensors is 40 mm. The sensor is set as an
ideal ultrasound detector with infinite bandwidth. The speed of sound outside
the skull in the simulation is set as 1500 m/s. These sensors recorded the PA
pressure of whole medium with 150 MHz sampling frequency.

Table 1: Speed, Density and Absorption of Distribution on Model According to Img-pixel-values Tissue | Speed | Density | Absorption | Img-pixel
---|---|---|---|---
| $[m/s]$ | ($Kg/m^{3}$) | $[dB/(MHz\cdot cm)]$ | -values
air | 343 | 1.2 | 0.0004 | [0]*
water | 1475 | 1000 | 0.05 | [0]*
midbrain | 1546.3 | 1000 | 0.6 | [21-39],[51-78]
white matter | 1552.5 | 1050 | 0.6 | [40-50]
grey matter | 1500 | 1100 | 0.6 | [81-220]
spinal fluid | 1475 | 1000 | 0.05 | [1-9]
scalp | 1540 | 1000 | 0.1 | [10-20]
skull | 3476 | 1979 | 2.7 | [221-255]

##  4 Results

###  4.1 TRPA Guided Focus within Distorted Medium

Since the speed of light is much faster than the speed of sound, ultrasound
can be regarded as static when the light passes through it, . so So here we
use the maximum signal pressure recorded at each point of the ROI to evaluate
the acoustic focusing ability. The PA signal was generated by the k-Wave
simulation in MATLAB, followed by transcranial transmission, and received by
64 acoustic sensors. We obtained the focus of the sound waves at the initial
position with and without TR. An absorber with size of 2 mm and absorption
coefficient of 5 $cm^{-1}$ is placed in the center of the ROI to generate a PA
signal. Experimental results are shown in Fig. 5. Our simulation environment
is carried out in a 3D configuration. In order to get the acoustic focus in
the space more intuitively, we intercepted three different planes of $xy$,
$xz$ and $xz$ for comparison. Using the focused transducer, it resulted in an
ultrasound focus at the intended location shown in Fig. 5(a). We found that
ultrasound is focused at the intended location, but there were many “hot
spots” outside the discs with a similar or even higher pressure gain than the
value inside focus. Local “hotspots” at the edge of the disc that extended out
into the surrounding tissue may lead to generate unwilling ‘guide star’ that
affects the accuracy of the tagging. In comparison, using TRPA guided focus
(Fig. 5(b)) resulted in an ultrasound focusing at the intended location with
much less “hotspots” on the disc outside.

Figure 5: Simulated pressure gain distribution in the cross-section plane of
xy, xz and yz. The black contours represent the shape of the skull. (a)
ultrasound direct focus. (b) TRPA guided focus. Figure 6: Simulated result of
mean intensity fluence images computed along the z-axis. (a) Without wavefront
shaping. (b) TRUE. (c)TRAP-TRUE.

###  4.2 TRPA-TRUE Optical Focus within Distorted Medium

Figure 6(a) shows the mean intensity projection of the 3D stack acquired
without wavefront shaping. We intercepted three cross-sections along y-axis to
more intuitively see the focus of light in space. It can be clearly seen from
Fig. 6(a) that the intensity of the flux decreases rapidly as the distance
from the scattering layer increases.

Figure 6(b) shows the mean intensity projection of the fluence image after
applying the wavefront shaping obtained from traditional TRUE method. Here we
can clearly see that TRUE greatly improves the light focusing effect compared
to that without wavefront shaping. However, around the middle of the z-axis
(i.e. position of ultrasonic ’guide star’), the focus of the light seems a bit
blurred with expanded focal size and one extra unwanted focal spot on top of
it. More severely, the highest intensity is around the surface of the tissue,
rather than at the focal spot It can also be seen that due to the low signal-
to-noise ratio, the intensity does not follow a monotonic change at these
depths.

Figure 6(c) shows the mean intensity projection of the fluence image after
applying wavefront shaping obtained from TRPA-TRUE method. We use the time-
reversed PA signal as an excitation of ultrasound modulation to tag the
propagating photons. In contrast to the traditional TRUE method, the fluence
heat map shows a much smaller focus area with higher intensity projection at
the target position.

##  5 Discussions And Conclusion

In this work, we demonstrated TRPA-TRUE by integrating photoacoustic,
ultrasonic time-reversal, and OPC, for the first time, achieving ultrafast
2D/3D optical focusing in both optically and acoustically heterogeneous medium
with one snapshot. Using the MC simulation and k-Wave toolbox, we demonstrated
the ability to enhance focus intensity with higher SNR through numerical brain
model, which has not been previously reported for any TRUE system. Although
time-reversed photoacoustic signals has shown the feasibility to provide
ultrasound focus for acousto-optic interaction in this work, for in vivo human
brain, the presence of the skull makes it challenging to generate intracranial
photoacoustic signals due to the strong light scattering and attenuation of
the skull, which is an open challenge to be solved.

##  6 Funding

This research was funded by Natural Science Foundation of Shanghai
(18ZR1425000), and National Natural Science Foundation of China (61805139).

##  7 Acknowledgments

The MRI data is downloaded from Neuromorphometrics, Inc[30] and the ultrasound
pressure data is downloaded on github from Yujia Huang. Thanks for their
support.

## References

  * [1] D. G. Grier, “A revolution in optical manipulation,” nature 424, 810–816 (2003).
  * [2] E. S. Boyden, F. Zhang, E. Bamberg, G. Nagel, and K. Deisseroth, “Millisecond-timescale, genetically targeted optical control of neural activity,” Nature neuroscience 8, 1263–1268 (2005).
  * [3] M. F. Yanik, H. Cinar, H. N. Cinar, A. D. Chisholm, Y. Jin, and A. Ben-Yakar, “Functional regeneration after laser axotomy,” Nature 432, 822–822 (2004).
  * [4] Y. M. Wang, B. Judkewitz, C. A. DiMarzio, and C. Yang, “Deep-tissue focal fluorescence imaging with digitally time-reversed ultrasound-encoded light,” Nature communications 3, 1–8 (2012).
  * [5] V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nature methods 7, 603–614 (2010).
  * [6] Y. Liu, C. Ma, Y. Shen, J. Shi, and L. V. Wang, “Focusing light inside dynamic scattering media with millisecond digital optical phase conjugation,” Optica 4, 280–288 (2017).
  * [7] J. Kubby, S. Gigan, and M. Cui, _Wavefront shaping for biomedical imaging_ (Cambridge University Press, 2019).
  * [8] A. Thendiyammal, G. Osnabrugge, T. Knop, and I. M. Vellekoop, “Model-based wavefront shaping microscopy,” arXiv pp. arXiv–2002 (2020).
  * [9] A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, “Controlling waves in space and time for imaging and focusing in complex media,” Nature photonics 6, 283–292 (2012).
  * [10] M. F. Yanik and S. Fan, “Time reversal of light with linear optics and modulators,” Physical review letters 93, 173903 (2004).
  * [11] Z. Yaqoob, D. Psaltis, M. S. Feld, and C. Yang, “Optical phase conjugation for turbidity suppression in biological samples,” Nature photonics 2, 110–115 (2008).
  * [12] I. M. Vellekoop and A. Mosk, “Focusing coherent light through opaque strongly scattering media,” Optics letters 32, 2309–2311 (2007).
  * [13] C. Yang, H. Lan, and F. Gao, “Deep learning enables robust and precise light focusing on treatment needs,” arXiv preprint arXiv:2008.06975 (2020).
  * [14] P. Lai, L. Wang, J. W. Tay, and L. V. Wang, “Photoacoustically guided wavefront shaping for enhanced optical focusing in scattering media,” Nature photonics 9, 126–132 (2015).
  * [15] R. Horstmeyer, H. Ruan, and C. Yang, “Guidestar-assisted wavefront-shaping methods for focusing light into biological tissue,” Nature photonics 9, 563–571 (2015).
  * [16] I. M. Vellekoop, E. Van Putten, A. Lagendijk, and A. Mosk, “Demixing light paths inside disordered metamaterials,” Optics express 16, 67–80 (2008).
  * [17] X. Xu, H. Liu, and L. V. Wang, “Time-reversed ultrasonically encoded optical focusing into scattering media,” Nature photonics 5, 154–157 (2011).
  * [18] K. Si, R. Fiolka, and M. Cui, “Breaking the spatial resolution barrier via iterative sound-light interaction in deep tissue microscopy,” Scientific reports 2, 748 (2012).
  * [19] Y. Suzuki, J. W. Tay, Q. Yang, and L. V. Wang, “Continuous scanning of a time-reversed ultrasonically encoded optical focus by reflection-mode digital phase conjugation,” Optics letters 39, 3441–3444 (2014).
  * [20] H. Ruan, M. Jang, B. Judkewitz, and C. Yang, “Iterative time-reversed ultrasonically encoded light focusing in backscattering mode,” Scientific reports 4, 7156 (2014).
  * [21] E. Bossy, K. Daoudi, A.-C. Boccara, M. Tanter, J.-F. Aubry, G. Montaldo, and M. Fink, “Time reversal of photoacoustic waves,” Applied physics letters 89, 184108 (2006).
  * [22] J. Zhang, H. Lan, Y. Ma, Y. Wang, F. Gao, and F. Gao, “Photoacoustic-guided transcranial hifu with combined time-reversal and genetic algorithm,” in _2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC),_ (IEEE, 2020), pp. 1915–1918.
  * [23] W. Leutz and G. Maret, “Ultrasonic modulation of multiply scattered light,” Physica B: Condensed Matter 204, 14–19 (1995).
  * [24] L. V. Wang, “Mechanisms of ultrasonic modulation of multiply scattered coherent light: an analytic model,” Physical Review Letters 87, 043903 (2001).
  * [25] G. Mahan, W. Engler, J. Tiemann, and E. Uzgiris, “Ultrasonic tagging of light: theory,” Proceedings of the National Academy of Sciences 95, 14015–14019 (1998).
  * [26] L. V. Wang, “Mechanisms of ultrasonic modulation of multiply scattered coherent light: a monte carlo model,” Optics letters 26, 1191–1193 (2001).
  * [27] Y. Huang, M. Cua, J. Brake, Y. Liu, and C. Yang, “Investigating ultrasound–light interaction in scattering media,” Journal of Biomedical Optics 25, 025002 (2020).
  * [28] A. C. Tam, “Applications of photoacoustic sensing techniques,” Reviews of Modern Physics 58, 381 (1986).
  * [29] A. R. Funke, J.-F. Aubry, M. Fink, A.-C. Boccara, and E. Bossy, “Photoacoustic guidance of high intensity focused ultrasound with selective optical contrasts and time-reversal,” Applied Physics Letters 94, 054102 (2009).
  * [30] R. Bakker, P. Tiesinga, and R. Kötter, “The scalable brain atlas: instant web-based access to public brain atlases and related content,” Neuroinformatics 13, 353–366 (2015).
  * [31] B. E. Treeby and B. T. Cox, “k-wave: Matlab toolbox for the simulation and reconstruction of photoacoustic wave fields,” Journal of biomedical optics 15, 021314 (2010).
  * [32] L. Wang, S. L. Jacques, and L. Zheng, “Mcml—monte carlo modeling of light transport in multi-layered tissues,” Computer methods and programs in biomedicine 47, 131–146 (1995).
  * [33] S. Powell and T. S. Leung, “Highly parallel monte-carlo simulations of the acousto-optic effect in heterogeneous turbid media,” Journal of biomedical optics 17, 045002 (2012).
  * [34] T. S. Leung and S. Powell, “Fast monte carlo simulations of ultrasound-modulated light using a graphics processing unit,” Journal of biomedical optics 15, 055007 (2010).

ref


# Analyzing and Enhancing Closed-loop Stability in Reactive Simulation

Wei-Jer Chang, Yeping Hu, Chenran Li, Wei Zhan, and Masayoshi Tomizuka  *This
work was supported by Hong Kong Centre Logistics Robotics.
W.J. Chang, Y. Hu, C. Li, W. Zhan and M. Tomizuka are with the Department of
Mechanical Engineering, University of California, Berkeley, CA 94720 USA
[weijer chang, yeping hu, chenran li, wzhan<EMAIL_ADDRESS>
###### Abstract

Simulation has played an important role in efficiently evaluating self-driving
vehicles in terms of scalability. Existing methods mostly rely on heuristic-
based simulation, where traffic participants follow certain human-encoded
rules that fail to generate complex human behaviors. Therefore, the reactive
simulation concept is proposed to bridge the human behavior gap between
simulation and real-world traffic scenarios by leveraging real-world data.
However, these reactive models can easily generate unreasonable behaviors
after a few steps of simulation, where we regard the model as losing its
stability. To the best of our knowledge, no work has explicitly discussed and
analyzed the stability of the reactive simulation framework. In this paper, we
aim to provide a thorough stability analysis of the reactive simulation and
propose a solution to enhance the stability. Specifically, we first propose a
new reactive simulation framework, where we discover that the smoothness and
consistency of the simulated state sequences are crucial factors to stability.
We then incorporate the kinematic vehicle model into the framework to improve
the closed-loop stability of the reactive simulation. Furthermore, along with
commonly-used metrics, several novel metrics are proposed in this paper to
better analyze the simulation performance.

##  I INTRODUCTION

Self-driving vehicles have the potential to build efficient and safe
transportation in society. To enable self-driving vehicles to safely and
efficiently navigate in human-involved environments, the design of interactive
algorithms for different modules (e.g., prediction, decision-making) is of
great importance. There are three general approaches to evaluate the overall
performance of self-driving vehicles: road-testing, heuristic-based
simulation, and log replay. Road testing is the most direct way to evaluate
the performance of autonomous vehicles in the real world. Nevertheless, the
testing process is expensive and time-consuming. In addition, it is hard to
directly evaluate the performance of autonomous vehicles under safety-critical
scenarios in the real world. For heuristic-based simulation, traffic
participants will follow specific human-encoded or traffic rules. Although
heuristic-based simulation can mimic basic driving behaviors, it is still
challenging to generate irregular and complex human behaviors. For log replay
approach, traffic participants will follow their recorded actions to imitate
real-world scenarios. However, traffic participants that are operated by the
log replay cannot react to the self-driving vehicle’s new actions that differ
from its recorded actions. For example, when a self-driving vehicle determines
to slow down (new action) instead of maintaining the same velocity (recorded
action) in the simulation, it might collide with any simulated vehicles
behind. To deal with the aforementioned problems, several recent works
proposed the concept of reactive simulation [1][2] to tackle the challenge of
bridging the human behavior gap between simulation and real-world scenarios.
The term ’reactive simulation’ means traffic participants in the simulation
environment will react to the self-driving vehicle similar to what human
drivers will do in the real world. Each traffic participant, or so-called
simulated agent, is simulated by reactive simulation, and the self-driving
vehicle is the vehicle we want to evaluate in the simulation environment. To
simulate realistic human behaviors, reactive simulation is usually implemented
by using a predictor that acts as a controller to control each simulated
agent’s actions. However, as we will show later in our experiments, if a
predictor is directly used in a closed-loop setting to consequently generate
simulated actions, simulated traffic participants will generate undesirable
behaviors. For example, the simulated traffic participants will easily drive
off-road or drive with jerky trajectories after few steps of iteration. We
regard such model performance as a lack of stability. To the best of our
knowledge, no work has explicitly discussed and analyzed the stability of the
reactive simulation frameworks. In this paper, an in-depth stability analysis
of the reactive simulation and a possible solution to improve stability are
provided. Our contributions are summarized as follows:

  1. 1.

Incorporate the kinematic vehicle model into the reactive simulation
framework.

  2. 2.

Provide a thorough stability analysis of the reactive simulation framework and
discuss the effects of the kinematic vehicle model.

  3. 3.

Propose several novel metrics to analyze the simulation performance better.

The remainder of this paper is structured as follows. The related works are
presented in section II. The closed-loop reactive simulation framework is
explained in Section III. In Section IV, experimental results and ablation
study are provided and discussed. Section V concludes the paper.

##  II Related works

###  II-A Scenario Generation

Scenario generation often uses generative models, such as Generative
Adversarial Network (GAN) and Variational Autoencoders (VAE), to simulate
driving scenarios in autonomous driving [3][4]. In addition, some controllable
generation methods [5][6] are proposed recently. [5] proposed a style
generative model that can generate diverse interactions with different safety
levels by altering style coefficients and SceGene [6] presented a bio-inspired
traffic scenario generation method in a controllable manner. However, scenario
generation only simulates the interactions based on historical states, and it
cannot simulate the reactive behaviors during the simulation.

###  II-B Reactive Behavior Simulation

There are two different reactive behavior simulation approaches: heuristic-
based and learning-based simulation. In the heuristic-based simulation, the
traffic participants are controlled by specific human encoded rules. A general
lane-changing car-following behavior model is proposed in [7]. The intelligent
driver model (IDM) [8] is a continuous microscopic car-following model that
has smooth braking and acceleration strategies. An overview of microscopic
traffic modeling and simulation techniques is provided in [9]. Recently, high-
fidelity simulators such as CARLA [10] and SUMO [11] have been widely used to
evaluate self-driving vehicle algorithms. Although the heuristic-based
approach can produce reasonable trajectories, it loses the diversity to
capture human behavior since human behavior is complex and hard to describe by
hand-crafted rules. For example, it is difficult to describe complex behavior
when human drivers negotiate during merging scenarios.

As a result, several learning-based reactive simulation works were proposed to
capture realistic multi-agent behaviors [1][2]. Both TrafficSim [1] and SimNet
[2] demonstrated the reactive simulation concept, which uses a neural network
as a predictor to control the traffic participants at each time step, and each
participant will ”react” to each other every time step. The main advantage of
the reactive simulation is that every agent’s new action will be affected by
other agents’ motion states to model the interaction. Usually, reactive
simulation is implemented by the ”closed-loop” concept: the future states are
continuously predicted based on the new states to simulate how agents will
interact with each other in the next few seconds. In addition, TrafficSim
demonstrated the usage of simulations as an effective data augmentation
approach for training a better motion planner, and SimNet utilized the
simulation model to evaluate an existing motion planner. The main difference
between scenario generation and reactive simulation is that scenario
generation will generate the driving scenarios directly based on historical
states (open-loop). On the other hand, in reactive simulation, the simulated
results will be based on the interaction between the self-driving vehicle and
the traffic participants at each time step (closed-loop).

Figure 1: Visualization of our proposed closed-loop simulation framework.
Based on map information and historical states, (a) the predictor backbone
will first predict agents’ control actions. Then, the control actions are fed
into (b) the kinematic layer to generate kinematically feasible trajectories.
At last, the smoothing layer will produce relatively smooth trajectories. In
each iteration, each agent will follow their simulated reactive actions, and
their motion states will be updated. The framework will iteratively predict
agents’ future trajectories based on the new states to simulate agents’ future
plans.

###  II-C Kinematics in Motion Prediction

Accurately predicting the motion and interaction of the surrounding vehicles
are necessary capabilities of autonomous vehicles. In recent years, machine
learning methods have greatly improved the capability of predicting the
traffic participants in different driving scenarios [12][13][14]. Nonetheless,
pure machine learning methods cannot guarantee the dynamic feasibility [15] of
trajectories. For example, the predicted vehicle may turn in a small radius
that does not satisfy the non-holonomic constraints of the vehicle. Therefore,
the concept of adding a kinematic layer to machine learning methods was
proposed [16][17]. [16] proposed a novel generative model structure to enable
neural networks to satisfy the kinematic constraints automatically, and [17]
combined a deep CNN model with the kinematic bicycle model for kinematic
feasibility guarantees. Recently, several works have incorporated planning
methods in neural networks to generate physically feasible trajectories [18]
[19]. [18] proposed a model-based generator to generate future trajectories
under explicit constraints and used a learning-based selector to capture
implicit interactions. Compared to adding a kinematic layer, these planning-
based methods can also include environmental constraints [18], but they may
lose the complexity of capturing short-term human behaviors.

##  III Closed-loop simulation framework

In this section, we present our closed-loop simulation framework, as shown in
Fig. 1. Closed-loop simulation means the simulator will iteratively simulate
traffic participants’ future plans based on the new simulated states. The
framework can be regarded as a controller that controls simulated agents’
motion. The main motivation behind the closed-loop concept is to simulate the
interaction between the traffic participants in every time step. Our framework
mainly consists of three components: a backbone predictor, a kinematic layer,
and a smoothing layer. The predictor will first predict agents’ control
actions based on map information and historical trajectories. Then, the
control actions are fed into the kinematic layer to generate feasible
trajectories. At last, the smoothing layer will produce relatively smooth
trajectories.

Figure 2: Illustration of the updating strategy of closed-loop simulation. At
time $t$, the model will predict the predicted agent’s states $S^{t:t+T}$.
Then, the closed-loop simulation framework will update the simulated reactive
actions $S^{t:t+T_{u}}$.

###  III-A Closed-loop Formulation

Based on map information $M$ and agents’ historical motion states
$S^{t-T_{obs}:t}$, we aim to simulate agents’ future states in the next
$T_{sim}$ time steps $S^{t:t+T_{sim}}$. Here, $T_{obs}$ denotes the time
length of observed historical states and $t$ is the current time step.
Specifically, at a given time step $t$,
$S^{t}={\left\\{s_{1}^{t},s_{2}^{t},...,s_{N}^{t}\right\\}}$ represents all
$N$ agents’ dynamic states, and the $k$-th actor’s state $s_{k}^{t}$ contains
its position $\left(x_{k}^{t},y_{k}^{t}\right)$, heading angle $\psi_{k}^{t}$,
velocity $v_{k}^{t}$, vehicle length $l_{k}$, and vehicle width $w_{k}$. The
updating strategy in closed-loop simulation is shown in Fig. 2. At each
iteration, the model will predict simulated agents’ future states $S^{t:t+T}$
in the next $T$ time steps. Then, the simulated agents will follow their
simulated reactive actions $S^{t:t+T_{u}}$ depending on the updated time
length $T_{u}$. Therefore, the new historical states
$S^{t^{\prime}-T_{obs}:t^{\prime}}$ at time $t^{\prime}$ in the next iteration
will be the concatenation of past historical states and the new states
$S^{t^{\prime}-T_{obs}:t^{\prime}}=\left[S^{t+T_{u}-T_{obs}:t},S^{t:t+T_{u}}\right]$.
The simulation framework can be served as an evaluation platform for motion
planners. In the simulation environment, the ego vehicle is controlled by a
motion planning algorithm, while other agents are closed-loop controlled by
the simulator [2]. The main advantage of this formulation is that we can
construct highly interactive and challenging scenarios to improve testing
efficiency. In addition, it can also enhance machine learning models’
prediction performance from a data augmentation perspective.

###  III-B Predictor

In this paper, we utilize Vectornet [20] as our backbone predictor in the
overall simulation framework, which is treated as a controller to simulate
each agent’s future behaviors. Instead of using raster-based images to
represent map information, Vectornet proposed vectorized representation to
encode both static HD map information and agents’ dynamic states as sequences
of vectors and used Graph Neural Network (GNN) to model interactions between
different entities. Note that the backbone predictor module can be replaced by
any other prediction model.

###  III-C Kinematic Layer

During the closed-loop simulation, we observed that if the predicted
trajectories had a small discrepancy between past trajectories, the model
would easily generate unreasonable behaviors, such as driving off-road or
driving with jerky trajectories. The detailed discussions can be found in
Section IV. Therefore, we use a kinematic layer [17] to generate kinematically
feasible trajectories by following the kinematic bicycle model relation from
the paper [15]:

$\displaystyle\dot{x}$ $\displaystyle=v\cos\left(\psi+\beta\right),$ (1)
$\displaystyle\dot{y}$ $\displaystyle=v\sin\left(\psi+\beta\right),$
$\displaystyle\dot{\psi}$
$\displaystyle=\frac{v}{l_{r}}\sin\left(\beta\right),$ $\displaystyle\dot{v}$
$\displaystyle=a,$ $\displaystyle\beta$
$\displaystyle={tan}^{-1}{\left(\frac{l_{r}}{l_{r}+l_{f}}\tan{\gamma}\right)}$

The relation of kinematic bicycle model parameters is shown in Fig. 1.
$\gamma$ is the steering angle, and $l_{r}$ and $l_{f}$ is the distance from
the center to the rear wheel and front wheel, respectively. $\beta$ is the
angle of current velocity $v$ of the center of mass with respect to the
vehicle’s longitudinal axis. $a$ denotes the acceleration of the vehicle. The
hyper-parameter $l_{r}$ in Eq. (1) denotes the distance from the rear wheel to
the center of the vehicle. If not provided in the data, $l_{r}$ can be
estimated by the historical states of the vehicle. In our settings, we use the
predictor backbone to estimate the ratio of $\frac{l_{r}}{l_{r}+l_{f}}$ and
calculate $l_{r}$ based on the vehicle’s length ${l_{r}+l_{f}}$ from data.

In order to output a sequence of actions, we utilize the aforementioned
equations to define a Recurrent Neural Network (RNN) unit in the kinematic
layer. Instead of directly outputting agents’ trajectories, we let the last
layer of the neural network output control parameters
$u_{t}=[a_{t},\beta_{t}]$. Note that we use $\beta_{t}$ instead of
$\gamma_{t}$ as the second input since $\beta_{t}$ and $\gamma_{t}$ have
static one-to-one relationships, and we noticed an improvement of the model
performance using $\beta_{t}$. By taking the last frame of historical states
as the initial states and the control parameters at every frame $t$, the RNN
unit will roll out each agent’s future trajectories and update its states.
This approach can guarantee that the model will output kinematically feasible
trajectories. In addition, the kinematic layer also incorporates the time
dependencies between each frame in the model. One of the main advantages of
the kinematic layer is that it can not only predict 2D positions
$\left(x^{t},y^{t}\right)$, but also predict velocity $v^{t}$ and heading
angle $\psi^{t}$ at the same time. Such characteristic enables direct
operations of the closed-loop simulation without estimating velocity $v^{t}$
and heading angle $\psi^{t}$ based on 2D positions. All of the parameters in
the kinematic layer are non-trainable, and we can regard the layer as
nonlinear functions that decode agents’ trajectories. The backbone predictor
combined with the kinematic layer can be viewed as predicting control
parameters in the kinematic bicycle model space instead of the x-y space.

The reason we want to incorporate the kinematic bicycle model into the
framework is to ensure the feasibility and smoothness of predicted
trajectories, and we will discuss the importance of smooth trajectories in the
closed-loop simulation. With the added kinematic layer, we can scale the
output range of the predictor to bound the acceleration and angle in
reasonable ranges. Besides using the kinematic bicycle model relation, we also
applied kinematic motion equations of particles in the kinematic layer with
$a_{x}$ and $a_{y}$ as the control parameters, which we called the simplified
model kinematic layer. Using kinematic motion equations of particles can be
regarded as simplifying the vehicle to a point mass. We then use the
accelerations $a_{x}$ and $a_{y}$ to roll out the agent’s trajectories with
the same RNN concept. This simplified kinematic layer seems to have the best
performance in terms of stability among three different settings: predictor
without kinematic layer, predictor with kinematic bicycle layer, and predictor
with simplified kinematic layer, and we will provide a detailed analysis in
the experiment section.

###  III-D Smoothing Layer: Weighted Average

To prevent sudden behavior switches during the closed-loop simulation process,
we propose a weighted average smoothing layer to eliminate trajectory’s
discontinuity when predicted states are added to the input during closed-loop
simulation. We denote the weighted average trajectory at time $t$ as
$\left[x^{t:t+T},y^{t:t+T}\right]_{w}^{T}$ and the new predicted trajectory at
time $t^{\prime}$ as
$\left[x^{t^{\prime}:t^{\prime}+T},y^{t^{\prime}:t^{\prime}+T}\right]_{p}^{T}$.
We use Eq. (2) to update the weighted average trajectory
$\left[x^{t:t+T},y^{t:t+T}\right]_{w}^{T}$. Then, the velocity and heading
angle are derived based on the weighted average trajectory. This method can
mitigate the discrepancy between the training and testing environment, which
is effective in ensuring the stability of the closed-loop simulation.

$\begin{bmatrix}x^{t^{\prime}:t^{\prime}+T}\\\
y^{t^{\prime}:t^{\prime}+T}\end{bmatrix}_{w}=(1-\alpha)\begin{bmatrix}x^{t^{\prime}:t^{\prime}+T}\\\
y^{t^{\prime}:t^{\prime}+T}\end{bmatrix}_{p}+\alpha\begin{bmatrix}x^{t:t+T}\\\
y^{t:t+T}\end{bmatrix}_{w}$ (2)

##  IV Experiments

In this section, we first showcase the failure cases by directly using a
predictor to predict 2D positions in the closed-loop simulation framework. We
showcase both visualization and quantitative analysis in the ablation study
for the kinematic layer and the smoothing layer in closed-loop simulation. The
sensitivity issue of the kinematic layer and the effect of the weighted
average smoothing layer are discussed.

###  IV-A Dataset

We train and test our model using the USA_Roundabout_FT data from the
INTERACTION dataset [21], which contains highly interactive driving cases. As
shown in Fig. 1, the roundabout has six branches, and each branch has two
directions (in and out), which creates more interactions in the scenario. In
this paper, we aim to test the closed-loop simulation approach and evaluating
the framework’s performance under this challenging scenario.

###  IV-B Implementation Details

For the predictor, we vectorized the HD map and past 1-second (10 frames)
agent trajectories and trained the VectorNet backbone to predict the next 3
seconds of the target agent’s trajectory. In the simulation, the simulated
agent has access to other agents’ states within an observable radius of 70
meters. For the simulation setup, we use 1-second ground truth data to
initialize the agent’s states and unroll the simulation for 5 seconds. We set
5-second as the simulation duration time since most of the interaction
scenarios covered in the dataset happen in this range. The updated frequency
is 10 Hz, which means we simulate 0.1 seconds in every iteration. At each
iteration, the simulated agent will take the control action at the first step,
and the new states will be the next input to the closed-loop simulation
framework.

To thoroughly analyze the closed-loop simulation framework, we set up three
different framework settings: predictor without kinematic layer, predictor
with kinematic bicycle layer, and predictor with simplified kinematic layer.
Each setting is composed of two sub-settings: with and without a weighted
average method. A total of six different framework settings are considered,
which details are shown below:

  * •

_xy_ : Directly use Vectornet as the predictor to simulate future
trajectories.

  * •

_xy weighted_ : Directly use Vectornet with a weighted average smoothing layer
to simulate future trajectories.

  * •

_kinematic_ : Use Vectornet to predict control parameters
$u_{t}=[a_{t},\beta_{t}]$, which will be fed into the kinematic layer with the
kinematic bicycle model equations to simulate future trajectories.

  * •

_kinematic weighted_ : Use Vectornet to predict control parameters
$u_{t}=[a_{t},\beta_{t}]$, which will be fed into the kinematic layer with the
kinematic bicycle model equations to simulate future trajectories with
weighted average smoothing.

  * •

_axay_ : Use Vectornet to predict control parameters $u_{t}=[a_{x},a_{y}]$,
which will be fed into the kinematic layer with the particle acceleration
equations to simulate future trajectories.

  * •

_axay weighted_ : Use Vectornet to predict control parameters
$u_{t}=[a_{x},a_{y}]$, which will be fed into the kinematic layer with the
particle acceleration equations to simulate future trajectories with weighted
average smoothing.

Note that to run the closed-loop simulation, we need to calculate the updated
motion state from the framework output. For _xy_ and _xy weighted_ , we use
Cubic spline [22] to interpolate the trajectories and use the interpolated
function to derive heading angle and velocity. Similarly, for _axay_ and _axay
weighted_ , we obtain the heading angle from Cubic spline interpolation, while
the updated velocity is directly obtained from the kinematic layer. To analyze
the model’s performance, we simplify the simulation process by only simulating
one car in each case and letting the other traffic participants follow their
recorded actions. For the smoothing layer, we choose $\alpha=0.2$.

###  IV-C Metrics

Since there are no standard evaluation metrics for closed-loop simulation, we
follow the metrics defined in [1][2] while proposing motion smoothness (MS)
and predicted trajectory difference (TD) to analyze the stability of the model
performance.

Collision Rate (CR): Collision rate reflects partial interaction of how the
simulated agent reacts to other traffic participants. We define collision rate
as the number of times the simulated agent collides out of the total number of
sampled scenarios. The simulated agent is considered to have a collision if
its bounding box overlaps with that of the other agents.

Simulation Realism: Average distance error (ADE) and final distance error
(FDE) between ground truth and simulated positions are reported to measure how
well the simulation can reconstruct real-world traffic scenarios.

Motion Smoothness (MS): During the closed-loop simulation, the smoothness is
strongly related to stability. Therefore, we use the absolute value of jerk,
the rate that the agent’s acceleration changes with respect to time to measure
the smoothness of simulated trajectories.

Trajectory Difference (TD): Trajectory difference is defined as the mean
square error between the predicted trajectories at a different time step.
Instead of the above closed-loop metrics, we also define open-loop metrics
(i.e., how the prediction changes at each iteration). Specifically, in our
setting, we predict 30 frames at each time step, so the difference is the
overlapped 29 frames between the two predictions. Trajectory difference
measures the consistency of the model when unrolling in closed-loop
simulation. Usually, an abrupt change in the ego vehicle’s behavior is
abnormal and will cause the predictor to generate unreasonable outputs.

Prediction metrics: We also provide the prediction metrics of each model,
including average distance error (ADE), final distance error (FDE), and Miss
Rate (MR). We consider a case as missed when the difference between the
predicted vehicle’s position and the ground truth position at the final frame
is larger than the given lateral or longitudinal threshold defined in [21]. MR
is the percentage of missed cases over all cases.

###  IV-D Analysis of Failure Cases

In this subsection, we demonstrate how directly using the predictor would
cause the model to fail easily during the closed-loop simulation. In fact,
when the simulated ego vehicle continuously follows its predicted
trajectories, the discontinuity of the predicted states will cause the model
to yield unreasonable trajectories.

When the simulated agent is entering the roundabout shown in Fig. 3(a), we can
observe that the ego vehicle starts to deviate at $t=1.2s$. According to the
velocity and heading angle plots, we can see a discontinuity starting from
$t=1.2s$. After that, the simulated agent starts to deviate from the road. In
another example shown in Fig. 3(b), the simulated agent slows down and stops
in the middle of the roundabout. Similarly, we can observe discontinuity in
velocity and heading angle profile.

Figure 3: Visualization of failure cases in _xy_ setting: (a) Entering
roundabout scenarios; (b) leaving roundabout scenarios. The simulated agent
loses its stability after a few steps. The velocity and angle profile show the
discontinuity of states in closed-loop simulation. For visualization purposes,
the dots are plotted every 0.3 seconds. Figure 4: Ablation study of kinematic
layer and smoothing layer. We input a 1-second historical states to simulate
the agent for 5 seconds. The simulated trajectories, velocities, and angle
profiles are shown in the figure. The figure shows that _axay weighted_
performs the best in all scenarios. We can observe that the heading angle
$\psi$ profile of _kinematic_ has oscillation in (a) and (c), which implies
that the kinematic bicycle layer is sensitive to the initial heading angle.

Based on our observation, we find that the stability of the closed-loop
simulation is strongly related to the consistency between the predicted and
past sequences. Once the ego vehicle takes a step that does not smoothly
connect with the past state sequences, it becomes difficult for the model to
recover from the resulted state deviation, which cause the model generates
either off-road or jerky trajectory. Therefore, how to maintain the
consistency between the simulated states and historical states becomes a
crucial problem when running closed-loop simulations.

###  IV-E Ablation Study

In this subsection, we compare six different framework settings: _xy_ , _xy
weighted_ , _kinematic_ , _kinematic weighted_ , _axay_ , and _axay weighted_.
The influences of adding the kinematic layer and the weighted average
smoothing layer to the model are discussed. Furthermore, differences between
kinematic bicycle layer _kinematic_ and simplified kinematic layer _axay_ are
also discussed.

One key observation we find is that when the velocity $v$ and heading angle
$\psi$ do not vary much during the process, the closed-loop simulation will
gradually deviate from the road instead of performing unreasonable behaviors
right after a few steps. In Fig. 4(a), _xy_ simulates the ego vehicle stopping
in the middle of the roundabout with jerky velocity and angle profile. On the
other hand, _kinematic_ , _kinematic weighted_ , and _axay weighted_ all
simulate reasonable trajectories with smooth state sequences. Similar
observations can be seen in both Fig. 4(b) and Fig. 4(c).

TABLE I:  Comparison between six different framework settings. The prediction
metrics reflect the prediction performance of each predictor (open-loop). The
closed-loop evaluation metrics reflect simulation performance, where each
simulation is simulated for 5 seconds. Note that for prediction metrics, it is
unnecessary to add a weighted average layer because the predictors are not
evaluated in closed-loop.

Model | Prediction Metrics | Closed-Loop Simulation Metrics
---|---|---
ADE (m) |  FDE (m) |  MR ($\%$) | ADE (m) |  Overall ADE (m) |  FDE (m) |  CR ($\%$) |  MS (m/$s^{3}$) |  TD (m)
0-1s | 1-2s | 2-3s | 3-4s | 4-5s
_xy_ | 0.57 (0.19) | 1.41 (0.4) | 20.6 (4.6) | 4.37 (0.59) | 10.32 (1.39) | 15.78 (2.19) | 21.74 (2.68) | 27.68 (2.46) | 15.56 (1.81) | 22.78 (1.87) | 32.0 (11.0) | 207.65 (26.22) | 37.51 (4.02)
_xy weighted_ | 3.41 (0.79) | 6.63 (1.58) | 10 (2.01) | 13.8 (2.36) | 17.89 (2.65) | 10.08 (1.83) | 14.82 (2.07) | 24.0 (6.1) | 73.42 (22.75) | 27.37 (3.69)
_kinematic_ | 0.55 (0.11) | 1.46 (0.29) | 19.3 (2.6) | 2.28 (0.25) | 5.43 (1.65) | 9.51 (3.26) | 14.03 (4.38) | 18.68 (4.85) | 9.69 (2.74) | 15.23 (4.09) | 22.0 (6.7) | 2.51 (0.76) | 11.52 (3.89)
_kinematic weighted_ | 2.27 (0.2) | 4.59 (0.52) | 7.49 (0.93) | 11.11 (1.46) | 15.07 (2.09) | 7.88 (0.98) | 12.87 (1.74) | 18.0 (2.4) | 8.61 (3.86) | 11.95 (2.7)
_axay_ | 0.65 (0.03) | 1.79 (0.09) | 21.3 (2.9) | 1.87 (0.06) | 5.08 (0.13) | 9.67 (0.26) | 15.44 (0.51) | 21.92 (0.84) | 10.76 (0.35) | 16.21 (0.61) | 27.0 (7.3) | 8.49 (3.87) | 74.04 (19.55)
_axay weighted_ | 1.67 (0.05) | 3.2 (0.12) | 5.37 (0.05) | 8.17 (0.14) | 11.35 (0.36) | 5.93 (0.06) | 9.87 (0.17) | 14.0 (0.9) | 11.1 (3.22) | 47.91 (3.47)

####  IV-E1 Sensitivity of the Kinematic Layer

Since the kinematic layer calculates future trajectories based on the velocity
$v$ and heading angle $\psi$, it is susceptible to the input initial states.
For models that incorporate the kinematic bicycle layer (i.e., _kinematic_ and
_kinematic weighted_), both the heading angle and velocity are crucial to the
predicted sequences. For example, in Fig. 4(a) and 4(c), we can see there is
an obvious oscillating of the heading angle of the simulated agent in
_kinematic_. We believe that the sensitivity issue of the kinematic layer may
cause the simulated trajectories to slightly deviate from the road, as shown
in Fig. 4(b) and 4(c). Similarly, _axay_ is sensitive to the initial
velocities $v_{x}$ and $v_{y}$. For example, in Fig. 4(b) and 4(c), when the
model _axay_ predicts a faster speed, the simulated agent keeps accelerating
and eventually deviate from the road.

####  IV-E2 Effect of the Weighted Average Method

The weighted average method prevents the simulated agent from making a sudden
change of trajectories and thus improves the stability of the model. The
effect can be seen in the improvement from _xy_ to _xy weighted_ setting. In
Fig. 4(b), _xy weighted_ successfully simulates the agent entering the branch.
In Fig. 4(c), the simulated agent in _xy_ suddenly drives toward the
divisional island with a jerky trajectory, while _xy weighted_ gradually
deviates from the road. In both cases, the velocity and heading angle profile
in _xy weighted_ do not have a dramatic oscillation compared to _xy_. As for
framework settings with the kinematic bicycle layer, the simulated agent under
the _kinematic_ setting performs similarly to that under the _kinematic
weighted_ setting.

The weighted average method also maintains stability in simplified kinematic
layer setting _axay_. From Fig. 4(a), the simulated agent in _axay_ will
either accelerate to drive off-road (Fig. 4(b) and Fig. 4(c)) or decelerate
and deviate from typical trajectories (Fig. 4(a)). We argue that the weighted
average helps the predicting state sequences maintain in a certain range that
the model can preserve its stability. This characteristic allows the closed-
loop simulation framework to have a tolerance for slightly deviated
predictions. However, this is not a perfect solution for closed-loop
simulation since the weighted average will introduce inertia to the system.
Thus, the simulated agent cannot react immediately to other vehicles.

Among all six different framework settings, _axay weighted_ performs the best
in terms of stability and reasonableness. We believe there are three main
reasons. First, using the kinematic layer will produce a smoother state
sequence to reduce the discrepancy between the historical and simulated
states. Second, in terms of sensitivity, _axay_ is only sensitive to
velocities, while _kinematic_ is sensitive to both velocities and heading
angles. We argue that the heading angle will have more influence on the
simulated trajectories, and this may be the reason why _axay weighted_ works
better than _kinematic_. At last, the weighted average method allows the
closed-loop simulation framework to have a tolerance for predictions that
deviate from normal trajectories.

###  IV-F Quantitative Analysis of Proposed Approach

Table 1 shows the comparison of six different framework settings in randomly
selected 100 scenarios. Although the prediction performance under setting
_axay_ is slightly worse than that under setting _xy_ and _kinematic_ , _axay
weighted_ outperforms other framework settings overall in closed-loop
simulation metrics. Also, according to the standard deviation of all these
metrics, _axay weighted_ has the most consistent performances. In terms of
motion smoothness, _xy_ has the largest MS metric with 207.65 m/$s^{3}$. This
large MS value reflects that directly using the predictor to simulate
trajectories will easily cause instability and produce jerky trajectories.
After adding a kinematic layer, we can observe that MS is greatly improved in
both _axay_ and _kinematic_ settings.

Compared with _xy_ , the overall performance of _xy weighted_ is greatly
improved, especially for the motion smoothness metric. Similarly, we can also
observe the improvement of _axay_ compared to _axay weighted_ under most
metrics. Note that the TD value for both _xy_ and _axay_ settings is dropped
in terms of mean and standard deviation after adding the weighted average
method. This implies that the weighted average improves the consistency of the
predictions at each timestamp and thus improves the closed-loop simulation
performance. On the other hand, for _kinematic_ to _kinematic weighted_
settings, both MS and TD do not have improvements, while simulation realism
(ADE and FDE) and collision rate (CR) are still improved. One potential reason
is that the weighted average smoothing layer improves the model’s tolerance
when the model produces a deviated state sequence. In addition, the ADE and
FDE of _kinematic_ and _kinematic weighted_ are larger than those of _axay
weighted_. We believe this result comes from the sensitivity issue of the
kinematic bicycle model: a slight angle deviation would cause the simulated
agent to accumulate its error during closed-loop simulation.

###  IV-G Further Discussion

Currently, there are still a few cases that will gradually drive off-road even
in our best framework setting, _axay weighted_ (shown in Fig. 5). Most of
these failure cases are under a situation that will easily cause model
confusion. For example, in Fig. 5(a), the simulated agent has two driving
options: one is to leave the roundabout and drive towards an exit branch, and
the other is to stay in the roundabout. If we use a deterministic neural
network model such as Vectornet to predict the agent’s trajectory, the model
will tend to pick the average of two plans in order to minimize the training
loss. Therefore, introducing multi-modality to the agents’ future plans will
resolve this issue. Another potential reason is that the state error of the
simulated vehicle starts to accumulate when we simulate the agent in closed-
loop. At each time step, the simulated agent takes actions based on slightly
deviated trajectories and eventually drives off-road. This is a common
distributional shift problem in imitation learning. In the future, how to deal
with distribution shift would be an important topic in closed-loop simulation
for long-horizon simulation. In general, it is common to assume a high-
performance predictor would also perform well under the closed-loop setting.
However, in this work, we discover that it is not the case, and both the
smoothness and consistency of the simulated state sequences are important to
the stability of closed-loop simulation.

Figure 5: Visualization of simulated agent’s trajectories, velocity, and angle
profile. This figure shows the distribution shift problem in the long horizon.
For visualization purposes, the dots are plotted every 0.3 seconds.

##  V CONCLUSION

In this paper, a novel simulation framework and a thorough stability analysis
of reactive simulation are presented. Our experiments show that directly using
a predictor in the reactive simulation will easily generate off-road or jerky
trajectories. We find out that during the closed-loop simulation, the
smoothness and consistency of the simulated state sequences are crucial to
stability. The kinematic bicycle layer can improve the stability but suffer
from the sensitivity issue of initial states. On the other hand, the
simplified kinematic layer performs the best in terms of simulation realism
and stability. Furthermore, we demonstrate how our weighted average method
improves the stability in reactive simulation. For future work, we will tackle
the long horizon distribution shift problem and consider the multi-modality of
the traffic scenarios.

## References

  * [1] S. Suo, S. Regalado, S. Casas, and R. Urtasun, “Trafficsim: Learning to simulate realistic multi-agent behaviors,” in _2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2021, pp. 10 395–10 404.
  * [2] L. Bergamini, Y. Ye, O. Scheel, L. Chen, C. Hu, L. Del Pero, B. Osiński, H. Grimmett, and P. Ondruska, “Simnet: Learning reactive self-driving simulations from real-world observations,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 5119–5125.
  * [3] A. Demetriou, H. Alfsvåg, S. Rahrovani, and M. H. Chehreghani, “A deep learning framework for generation and analysis of driving scenario trajectories,” _CoRR_ , vol. abs/2007.14524, 2020. [Online]. Available: https://arxiv.org/abs/2007.14524
  * [4] W. Ding, M. Xu, and D. Zhao, “Cmts: A conditional multiple trajectory synthesizer for generating safety-critical driving scenarios,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 4314–4321.
  * [5] Z. Yin, L. Sun, L. Sun, M. Tomizuka, and W. Zhan, “Diverse critical interaction generation for planning and planner evaluation,” _CoRR_ , vol. abs/2103.00906, 2021. [Online]. Available: https://arxiv.org/abs/2103.00906
  * [6] A. Li, S. Chen, L. Sun, N. Zheng, M. Tomizuka, and W. Zhan, “Scegene: Bio-inspired traffic scenario generation for autonomous driving testing,” 2021, pp. 1–16.
  * [7] P. Gipps, “A behavioural car-following model for computer simulation,” _Transportation Research Part B: Methodological_ , vol. 15, no. 2, pp. 105–111, April 1981. [Online]. Available: https://ideas.repec.org/a/eee/transb/v15y1981i2p105-111.html
  * [8] M. Treiber, A. Hennecke, and D. Helbing, “Congested traffic states in empirical observations and microscopic simulations,” _Phys. Rev. E_ , vol. 62, pp. 1805–1824, Aug 2000. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevE.62.1805
  * [9] A. Kesting, M. Treiber, and D. Helbing, “Agents for traffic simulation,” in _Multi-Agent Systems_ , 2009.
  * [10] A. Dosovitskiy, G. Ros, F. Codevilla, A. M. López, and V. Koltun, “CARLA: an open urban driving simulator,” vol. abs/1711.03938, 2017. [Online]. Available: http://arxiv.org/abs/1711.03938
  * [11] P. A. Lopez, M. Behrisch, L. Bieker-Walz, J. Erdmann, Y.-P. Flötteröd, R. Hilbrich, L. Lücken, J. Rummel, P. Wagner, and E. Wiessner, “Microscopic traffic simulation using sumo,” in _2018 21st International Conference on Intelligent Transportation Systems (ITSC)_ , 2018, pp. 2575–2582.
  * [12] Y. Hu, W. Zhan, and M. Tomizuka, “Probabilistic prediction of vehicle semantic intention and motion,” in _2018 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2018, pp. 307–313.
  * [13] Y. Hu, X. Jia, M. Tomizuka, and W. Zhan, “Causal-based time series domain generalization for vehicle intention prediction,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 7806–7813.
  * [14] Y. Hu, W. Zhan, and M. Tomizuka, “Scenario-transferable semantic graph reasoning for interaction-aware probabilistic prediction,” _arXiv preprint arXiv:2004.03053_ , 2020.
  * [15] J. Kong, M. Pfeiffer, G. Schildbach, and F. Borrelli, “Kinematic and dynamic vehicle models for autonomous driving control design,” in _2015 IEEE Intelligent Vehicles Symposium (IV)_ , 2015, pp. 1094–1099.
  * [16] H. Ma, J. Li, W. Zhan, and M. Tomizuka, “Wasserstein generative learning with kinematic constraints for probabilistic interactive driving behavior prediction,” _2019 IEEE Intelligent Vehicles Symposium (IV)_ , pp. 2477–2483, 2019.
  * [17] H. Cui, T. Nguyen, F. Chou, T. Lin, J. Schneider, D. Bradley, and N. Djuric, “Deep kinematic models for physically realistic prediction of vehicle trajectories,” _CoRR_ , vol. abs/1908.00219, 2019. [Online]. Available: http://arxiv.org/abs/1908.00219
  * [18] H. Song, D. Luan, W. Ding, M. Y. Wang, and Q. Chen, “Learning to predict vehicle trajectories with model-based planning,” in _Proceedings of the 5th Conference on Robot Learning_ , ser. Proceedings of Machine Learning Research, A. Faust, D. Hsu, and G. Neumann, Eds., vol. 164. PMLR, 08–11 Nov 2022, pp. 1035–1045. [Online]. Available: https://proceedings.mlr.press/v164/song22a.html
  * [19] H. Girase, J. Hoang, S. Yalamanchi, and M. Marchetti-Bowick, “Physically feasible vehicle trajectory prediction,” _CoRR_ , vol. abs/2104.14679, 2021\. [Online]. Available: https://arxiv.org/abs/2104.14679
  * [20] J. Gao, C. Sun, H. Zhao, Y. Shen, D. Anguelov, C. Li, and C. Schmid, “Vectornet: Encoding HD maps and agent dynamics from vectorized representation,” _CoRR_ , vol. abs/2005.04259, 2020. [Online]. Available: https://arxiv.org/abs/2005.04259
  * [21] W. Zhan, L. Sun, D. Wang, H. Shi, A. Clausse, M. Naumann, J. Kümmerle, H. Königshof, C. Stiller, A. de La Fortelle, and M. Tomizuka, “INTERACTION Dataset: An INTERnational, Adversarial and Cooperative moTION Dataset in Interactive Driving Scenarios with Semantic Maps,” 2019. [Online]. Available: https://github.com/interaction-dataset/INTERPRET˙challenge˙single-agent
  * [22] S. Dyer and J. Dyer, “Cubic-spline interpolation. 1,” _IEEE Instrumentation Measurement Magazine_ , vol. 4, no. 1, pp. 44–46, 2001.