[2]Adeline Paiement

1]Université de Lille, Centrale Lille, Lille, France

2]Université de Toulon, Aix Marseille Univ, CNRS, LIS, Marseille, France

3]Observatoire de Paris/PSL, LESIA, CNRS, Paris, France

# Removing cloud shadows from ground-based solar imagery

Amal Chaoui<EMAIL_ADDRESS>Jay Paul Morgan
<EMAIL_ADDRESS><EMAIL_ADDRESS>Jean Aboudarham
<EMAIL_ADDRESS>[  [  [

###### Abstract

The study and prediction of space weather entails the analysis of solar images
showing structures of the Sun’s atmosphere. When imaged from the Earth’s
ground, images may be polluted by terrestrial clouds which hinder the
detection of solar structures. We propose a new method to remove cloud
shadows, based on a U-Net architecture, and compare classical supervision with
conditional GAN. We evaluate our method on two different imaging modalities,
using both real images and a new dataset of synthetic clouds. Quantitative
assessments are obtained through image quality indices (RMSE, PSNR, SSIM, and
FID). We demonstrate improved results with regards to the traditional cloud
removal technique and a sparse coding baseline, on different cloud types and
textures.

###### keywords:

Image cleaning, solar imaging, deep learning, U-Net, C-GAN

## Acknowledgments

This work was funded by ANR grant N° ANR-20-CE23-0014-01. One of the authors
(JA) thanks CNES (Centre National d’Etudes Spatiale, the French space agency)
for its support in this project.

##  1 Introduction and previous works

Figure 1: Overview of cloud removal from solar imagery. A cloud mask is
predicted and subtracted from the original image, to prepare for a subsequent
detection of solar structures (here ARs in red boxes). Our proposed method
only deals with the cloud removal stage.

This paper addresses the problem of cleaning ground-based solar imagery from
cloud shadow contaminants, as illustrated in Fig. 1, in support to forthcoming
analyses with solar physics and space weather applications. Monitoring the
activity of the Sun is a task of growing importance in the frame of space
weather research. The main sources of severe space weather disturbances at
Earth are energetic outbursts from the Sun, e.g. solar flares and coronal mass
ejections. They are strongly associated with solar structures of the
chromosphere (e.g. active regions (AR) and filaments, as illustrated in Fig.
5), visible from the Earth ground in the light emitted by ionised calcium (Ca-
II, at a wavelength of 3933.7Å) and the H-$\alpha$ line of neutral hydrogen
(6562.8Å). Despite their affordability, images taken from the ground present a
particular difficulty: they are prone to be tainted by cloud shadows when the
sky is not clear. These shadows are rarely uniform and may be mistaken for
solar structures (e.g. filaments). Thus, it is substantial to remove them
before automatic solar structure detection and monitoring.

Very few previous works addressed this problem. [1, 2] designed denoising
methods for H-$\alpha$ images. Feng et al. [1] proposed a metric for assessing
the level of cloud coverage in an image by considering the homogeneity of
image intensity between different quadrants of the solar disk. Additionally,
they proposed a denoising method that first estimates the pixel-wise
transmittance of clouds $\hat{I}_{t}$ from a cloud-free temporal neighbour
using a low-pass morphological filter to smooth out solar structures (i.e.
filaments). Second, it uses it to normalise the original image
$I_{\text{ini}}$:

$\hat{I}_{\text{clean}}=I_{\text{ini}}/\hat{I}_{t}$ (1)

In [2], Fuller and Aboudarham perform two stages of median filtering to
identify large scale intensity variations that correspond to clouds, while
thresholding out the solar structures. In both methods, the use of low-pass
filtering to distinguish clouds from solar structures is a possible limitation
when their spatial scales are similar in at least one direction.

In other application domains, some works removed large contaminants, such as
clouds and shadows, from remote sensing or natural images. Due to clouds
transparency, the task of removing clouds in Sun observations resembles to a
degree to these tasks. While there is a large literature on denoising methods
e.g. [3], we focus on the removal of large and semi-transparent contaminants,
as this is the closest to our task. Methods that perform inpainting over
opaque contaminants, such as [4] for thick clouds, are not considered to avoid
introducing artefacts in the study of solar structures (images with thick
clouds are therefore not usable). We also discard methods that exploit the
physics of image formation for natural images, e.g. [5, 6, 7, 8], as solar
imaging follows a different physics. Thin cloud and shadow removal methods are
good candidates for our task of removing semi-transparent clouds while
preserving the background solar structures.

Some shadow suppression methods e.g. [9] rely upon the sparse approximation
theory. A local dictionary of atoms is learnt from image patches as a basis
for sparse coding. This allows decomposing the image into reflectance
(invariant to shading), geometry (fine structural details), and shadow
components, the latter being excluded during image reconstruction. The use of
local dictionaries, rather than generic ones (e.g. wavelets, curvelets), helps
in learning atoms that adapt to each single image. However, this approach
cannot gather further insights from multiple examples.

This limitation is addressed by (deep) learning-based methods. [10, 11] use
encoder-decoder deep neural networks (DNN) to remove shadows by learning a
mapping between the initial images $I_{ini}$ and their binary mask $\hat{M}$
of shadow ratio, which is reminiscent of the transmittance image of [1]:

$\hat{I}_{\text{clean}}=I_{\text{ini}}/\hat{M}$ (2)

[12] and pix2pix [13] in [4] rely on skip connections (U-Net-style) to
preserve spatial resolution, and a conditional generative adversarial network
(C-GAN) for more efficient training. They obtain state-of-the-art (SoTA)
results in shadow and thin cloud removal. In [12], an additional C-GAN outputs
the denoised image from predicted shadow mask, while pix2pix of [4] outputs
the cleaned image directly.

We explore the applicability of these denoising methods for the new
application of removing cloud shadows from solar imagery. We evaluate sparse
representation-based and deep learning-based approaches, and compare against
[1, 2] (Section 4). Inspired by [12, 4], we use a combination of U-Net and
C-GAN (Section 3) to obtain the new SoTA for this task.

The methods are comparatively assessed on two imaging modalities. To support
these experiments, we introduce two new datasets of Sun images with artificial
and real clouds.

##  2 Data

Our images are from Paris-Meudon (PM) observatory111The Paris-Meudon
observatory is part of Observatoire de Paris-PSL:
https://www.observatoiredeparis.psl.eu/-observatoire-de-paris-.html?lang=en.
They are $1024\times 1024$ px full-disk spectroheliograms. Care has been taken
in their collection to ensure a good representativity. We used images of two
modalities: Ca-II (for studying ARs), and H-$\alpha$ (for filaments), as
highlighted in Fig. 5. These are major modalities in ground-based solar
imaging. Observations span different periods of the solar cycle in order to
reflect variations in the Sun’s activity and number / size of the solar
structures that may be observed at the surface of the Sun. 96 Ca-II and 111
H-$\alpha$ images are sampled between May 2002 and June 2017. Sampling uses
the metric of cloud coverage level proposed by [1] to select the images that
are tainted with cloud shadows.

PM images are taken in principle during clear sky days. However, even during
those days, small thin clouds emerging in the horizon cannot be prevented and
result in light cloud stains. Cloud shadows can come with different textures.
In PM observatory, the Sun disk is scanned by a thin moving slit over 20 to
100 sec. [14]. In case of fast enough cloud motion, due to the slit movement
alongside the motion of clouds, cloud spots may be shifted along the scanning
process, creating straight lines on the image (denoted as “streaked” texture).
Slower cloud motions result in the “fluffy” texture that is customary of many
terrestrial clouds. It is also a common texture in solar observations from
observatories that use 2D cameras rather than a spectrograph. The “streaked”
texture may provide more challenges to denoising methods, due to having
different spatial scales in different directions. Fig. 5 illustrates samples
of streaked and fluffy textures with different cloud cover thicknesses and
streak orientations.

Figure 2: Examples of real (top) and synthetic (bottom) cloud contamination in
Ca-II images.

The cloudy images don’t have cloud-free counterparts for training DNNs. Thus,
we introduce a dataset of real (initially cloud-free) images with added
synthetic cloud shadows. Using the cloud coverage level metric of [1], we
selected reasonably cloud-free images in the same time period as above. A
basic cloud texture was obtained from [15]. For each cloud-free image, a new
cloud shadow was composed by summing a random number (in range [2, 3]) of
duplicates of the basic texture with horizontal and vertical random resizing
between 50-100% of the original size, random flips (around one or two axes),
and random transparencies (in range [0.1, 0.4]). It was then overlapped to the
cloud-free image. This procedure provides realistic images of fluffy and semi-
transparent cloud shadow contamination. They are diverse and representative of
possible cloud shapes and opacity levels. The Ca-II and H-$\alpha$ synthetic
datasets comprise respectively 319 and 367 pairs of shadow/shadow-free images,
split into 179/44/96 and 205/51/111 training/validation/testing pairs. The
combination of our real and artificial cloud datasets provides a variety of
streaked and fluffy cloud textures. They are publicly
available222https://doi.org/10.5281/zenodo.7684200, along with the pretrained
cloud-removal DNNs333The source code and model weights are available via
https://gitlab.lis-lab.fr/presage/cloud-removal-from-solar-imagery..

Images have their intensity normalised in range $[0,1]$, centred solar disks
with normalised radii, and background pixels are set to zero outside of the
solar disk to remove any distracting features from background stars and
imaging artefacts.

##  3 Methodology

We adopt the deep learning approach due to its SoTA results in large
contaminants removal. As in [12, 13], we use skip connections to preserve
spatial resolution and details444We use the U-Net implementation from
https://github.com/mateuszbuda/brain-segmentation-pytorch.. Inspired by the
results of previous methods, we compare the classical fully supervised and the
C-GAN training strategies.

We experiment with our DNNs predicting either the cleaned image directly, or a
cloud shadow mask, to be used as in [1, 10, 11]. To avoid division by zero and
improve stability, we compute the clean image as:

$\hat{I}_{\text{clean}}=\frac{I_{\text{ini}}}{\hat{M}+\epsilon}$ (3)

We empirically verified that the exact value of $\epsilon$ has little impact
on the results, and we set it to $1e^{-5}$ in our experiments. Regressing a
shadow mask instead of the cleaned image directly is motivated by the physics
of solar imaging, as in [1]. In addition, estimating a transmittance or shadow
ratio, rather than the pixel’s absolute intensity, reduces the risk to create
artefacts that would hinder the forthcoming study of solar structures.
Finally, estimating a shadow mask should be a less complex task for the DNN,
which may be learnt more efficiently. Indeed, reproducing the exact patterns
of the more complex and diverse solar structures in the output image would be
an unnecessary additional task for the DNN.

With Eq. 3, the loss gradient varies as $I_{\text{ini}}/\hat{M}^{2}$, meaning
that the darkest pixels of the images, where cloud thickness is high, see
higher loss gradients and require a finer adjustment of their mask value. In
an attempt to improve stability during training through more consistent loss
gradient magnitudes, we also experiment with an alternative implementation,
i.e. adding in the intensity that was removed by clouds. This is similar to
predicting a residual mask:

$\hat{I}_{\text{clean}}=I_{\text{ini}}+\hat{M}$ (4)

We compare two training strategies: a fully supervised DNN (denoted as FS) and
a C-GAN following the pix2pix setup. For ease of comparison, the FS DNN and
C-GAN generator have the same U-Net-style architecture${}^{\ref{foot:unet}}$:
down-sampling and up-sampling blocks along with skip connections, each block
as defined in [16]. We experiment with a varying number of blocks to determine
the optimal complexity and receptive field needed to account for both large
scale contaminants and smaller scale solar structures. Due to memory
constraints, we go up to 6 down/up-sampling blocks, and find it to be the best
performing architecture (see supplementary materials). In the future, and on a
different hardware, more blocks could be tried.

We use a linear activation in the last block when predicting
$\hat{I}_{\text{clean}}$ directly or $\hat{M}$ as in Eq. 4. For Eq. 3, we find
that imposing the output to be in range $[0,1]$ is necessary to avoid
instability issues, so we use $0.5\times\tanh+0.5$ as activation.

The discriminator in the C-GAN is the same as in pix2pix [13]: a PatchGAN
network that considers image patches (small receptive fields) to capture high-
frequencies. We experimentally verified that the recommended receptive field
of $70\times 70$ is also optimal in our application. The discriminator is
applied to the cleaned image rather than the predicted shadow mask to assess
the preservation of solar structures.

For the C-GAN, we use the loss function suggested in [13]. The loss function
of the FS network and of the C-GAN generator uses the pixel-wise $L_{1}$-norm
since it promotes sparsity and encourages less blurring as noted in [13]. The
pixel-level loss is computed between the target $I_{t}$ and cleaned image
$\hat{I}_{\text{clean}}$:

$\mathcal{L}_{L_{1}}(I_{t},\hat{I}_{\text{clean}})=\frac{1}{N}\sum_{i=1}^{N}|I_{t}^{i}-\hat{I}_{\text{clean}}^{i}|$
(5)

For both DNNs, we use the default U-Net${}^{\ref{foot:unet}}$ training scheme
and batch size of 3. We train for 100 epochs and test on the epoch with the
best validation loss (FS DNN) or generator’s training loss (C-GAN, following
the usual practice for GANs). During testing of the C-GAN, the generator is
run with dropout, and batch normalisation considers the statistics of the
batch, as in [13].

##  4 Results

###  4.1 Choice of DNN output

We compare the three possible outputs of our DNNs, i.e. direct cleaned image
and shadow masks to be used as in Eqs. 3 (transmittance) and 4 (residual).
Table 1 reports the metrics on the Ca-II synthetic dataset and Fig. 4 provides
some illustration examples. For both DNNs, Eq. 3 is less stable than Eq. 4
during training, with poorer metrics and more noisy loss curves as shown in
Fig. 3. This more unstable training results in a poor convergence. Poorer
results are observed in the sensitive cases of strong cloud coverage (i.e. low
values in $I_{\text{ini}}$) where more subtle adjustments are needed in the
mask’s values to correct for the clouds. Predicting a residual mask $\hat{M}$
for Eq. 4 performs slightly better than the direct prediction of
$\hat{I}_{\text{clean}}$, which, in rare cases, created artefacts that may be
confused with solar structures. The better performance of Eq. 4 is in line
with other works on residual prediction. Following these observations, we
choose Eq. 4 to be used in the rest of this paper.

Table 1: Comparison of the three DNN outputs on the Ca-II synthesised dataset. Results are in format mean(std) over 10 runs. SSIM and RMSE are as $e-2$. DNN | Output | PSNR $\uparrow$ | SSIM $\uparrow$ | RMSE $\downarrow$ | FID $\downarrow$
---|---|---|---|---|---
Fully sup. | $\hat{I}_{\text{clean}}$ | 30.2(0.4) | 98.2(0.5) | 3.9(0.1) | 17.2(5.0)
$\hat{M}$ Eq. 3 | 28.9(0.2) | 98.8(0.1) | 4.5(0.1) | 20.2(2.4)
$\hat{M}$ Eq. 4 | 30.6(0.3) | 98.9(0.1) | 3.7(0.1) | 16.8(5.3)
C-GAN | $\hat{I}_{\text{clean}}$ | 29.4(0.6) | 97.7(0.4) | 4.1(0.2) | 22.0(10.2)
$\hat{M}$ Eq. 3 | 27.3(0.6) | 97.6(0.8) | 5.1(0.2) | 49.0(10.8)
$\hat{M}$ Eq. 4 | 30.0(0.4) | 98.8(0.2) | 3.9(0.1) | 18.4(6.3)
Figure 3: Loss curves when training the proposed models using Eqs. 3 and 4.
Figure 4: Comparative illustration of the three DNN outputs on some Ca-II
images with synthetic clouds. We display the results of a random run of each
network, since all runs produced roughly the same results.

###  4.2 Comparative evaluation of the methods

We compare the proposed DNNs against the SoTA [1, 2] and the sparse
approximation method of [9] as non-learning baselines. [1] requires a clean
temporal neighbour image (within 24 hours), which may not be always available.
We consider these cases a failure and report their numbers. The first
reflectance step of [9] is skipped, as reflection plays no role in solar
imaging. We use the Peak Signal-to-Noise Ratio (PSNR), Structured Similarity
Index Measure (SSIM) and Root Mean Squared Error (RMSE) metrics to assess the
quality of the predicted shadow-free images. In addition, the Fréchet
Inception Distance (FID) metric [17] is used to compare the distribution of
cleaned images against that of targets. For [1], the metrics are computed over
non-failure cases. For the DNNs, they are averaged over 10 runs.

Table 2: Performance on synthesised datasets. Results are in format mean(std) over 10 runs for the DNNs. SSIM and RMSE are as $e-2$. Data | Method | PSNR $\uparrow$ | SSIM $\uparrow$ | RMSE $\downarrow$ | FID $\downarrow$
---|---|---|---|---|---
Ca-II |  [1] (15 failures) | 21.90 | 94.8 | 8.6 | 104.8
[2] | 26.0 | 98.1 | 5.4 | 70.3
[9] | 23.2 | 92.7 | 7.3 | 45.5
Fully sup. | 30.6(0.3) | 98.9(0.1) | 3.7(0.1) | 16.8(5.3)
C-GAN | 30.0(0.4) | 98.8(0.2) | 3.9(0.1) | 18.4(6.3)
H-$\alpha$ |  [1] (5 failures) | 14.8 | 91.4 | 19.6 | 118.3
[2] | 23.3 | 98.4 | 7.0 | 134.7
[9] | 21.0 | 96.0 | 9.2 | 110.9
Fully sup. | 28.6(0.3) | 98.9(0.1) | 4.5(0.2) | 33.0(4.9)
C-GAN | 28.3(0.5) | 98.9(0.2) | 4.6(0.2) | 37.5(3.5)

Table 2 presents results on the synthetic cloud datasets, and Fig. 5
illustrates how the algorithms perform on real/synthetic clouds. More
illustrations are provided in the supplementary materials. [1] obtains the
worst results on both imaging modalities. It struggles to preserve solar
structures when their spatial scale is similar to the scale of clouds in at
least one direction, as often happens in the case of streaked textures. It is
also less robust to strong cloud thickness. [2] also struggles with the
spatial scale of streaked clouds. [9] fails to recover the underlying solar
structures and results in poor contrast.

The proposed DNNs outperform non-learning methods on all datasets. On
synthesised images, the DNN’s respective metrics are not significantly
different. Upon visual inspection of artificial and real images, FS produces a
slightly more uniform contrast across the solar disk, while the C-GAN’s
discriminator, focusing on small patches, allows more diversity of contrasts.
However, FS tends to over-compensate strong clouds resulting in a lower
contrast there, while C-GAN benefits from enforcing a patch-wise visually
pleasant output.

Figure 5: Examples of removal of real and synthetic clouds from Ca-II and
H-$\alpha$ images by compared algorithms. Images with real clouds (Real Ca-II,
Real H-$\alpha$, Ca-II AR, H-$\alpha$ Filament) do not have a clean target
counterpart, as explained in Section 2. Therefore, for these columns, the
target row is left blank.

It is worth noting that in our real-cloud dataset, most clouds are streaked,
while the DNNs are trained on the fluffy synthetic textures. Both DNNs
generalised well to this new cloud texture. While, in some cases, faint
streaks remain from the cleaned images, this happens for the strongest clouds
and generally doesn’t prevent the restoration of underlying solar structures.
C-GAN performs slightly better than FS in this particular case of strong
streaked clouds. We expect that in the future, training the DNNs with examples
of streaked clouds may help in handling the strongest streaked clouds.

Both DNNs provide a detailed restoration of ARs, with the aforementioned
contrast differences. Filaments suffer from a lower contrast, likely due to
their dark appearance, making them harder to distinguish from clouds. This is
confirmed by the thicker and therefore blurrier ones being the most affected
(erased) by this contrast reduction. C-GAN preserves slightly better these
thick filaments.

##  5 Conclusion

We investigated deep learning methods for removing large cloud shadow
contaminants from ground-based solar imagery. We adopted the U-Net-style
architecture for more detailed outputs. We compared three possible inference
methods (i.e. direct and through a transmittance or residual mask), and two
training setups: fully-supervised and C-GAN. The DNNs obtained better results
than the domain’s SoTA and a non-leaning method based on sparse approximation,
with better restoration of the underlying solar structures, and are the new
SoTA. The two training setups have different strengths and weaknesses: full
supervision produces more homogeneous solar disks that are more physically
plausible, and C-GAN better handles strong clouds and thick filament
restoration. Thus, there is no obvious “best setup”, and the choice of
training setup depends on the intended usage.

The DNNs have been trained on fluffy cloud textures only. While they
demonstrate a good generalisation to streaked cloud textures, further
improvement might be obtained from training on streaked clouds in the future.
Future work may also further focus on filament restoration, and evaluate the
effect of cloud removal on a subsequent image analysis e.g. the automatic
detection of solar structures.

Supplementary information

Supplementary materials are provided with additional comparative results.

Acknowledgments

This work was funded by ANR grant N° ANR-20-CE23-0014-01. One of the authors
(JA) thanks CNES (Centre National d’Etudes Spatiale, the French space agency)
for its support in this project.

## Declarations

The data used in this study is available at
https://doi.org/10.5281/zenodo.7684200. The code is available at
https://gitlab.lis-lab.fr/presage/cloud-removal-from-solar-imagery.

Adeline Paiement and Jean Aboudarham contributed to the study conception and
design. Data preparation was performed by Jay Paul Morgan, experiments were
performed by Amal Chaoui and Jay Paul Morgan, analysis was performed by Amal
Chaoui, Jay Paul Morgan and Adeline Paiement. All authors contributed to the
manuscript.

## References

  *   * Feng et al. [2014] Feng, S., Lin, J., Yang, Y., Zhu, H., Wang, F., Ji, K.: Automated Detecting and Removing Cloud Shadows in Full-Disk Solar Images. New Astronomy 32, 24–30 (2014)
  * Fuller and Aboudarham [2004] Fuller, N., Aboudarham, J.: Automatic detection of solar filaments versus manual digitization. In: Knowledge-Based and Intelligent Information and Engineering Systems (KES), pp. 467–475 (2004)
  * Chen et al. [2020] Chen, Z., Zeng, Z., Shen, H., Zheng, X., Dai, P., Ouyang, P.: DN-GAN: Denoising Generative Adversarial Networks For Speckle Noise Reduction in Optical Coherence Tomography Images. Biomedical Signal Processing and Control 55, 101632 (2020)
  * Zheng et al. [2021] Zheng, J., Liu, X.-Y., Wang, X.: Single Image Cloud Removal Using U-Net and Generative Adversarial Networks. IEEE Transactions on Geoscience and Remote Sensing 59(8), 6371–6385 (2021)
  * Finlayson et al. [2006] Finlayson, G.D., Hordley, S.D., Lu, C., Drew, M.S.: On the removal of shadows from images. IEEE Trans. Pattern Anal. Mach. Intell. 28, 59–68 (2006)
  * Yang et al. [2012] Yang, Q., Tan, K.H., Ahuja, N.: Shadow removal using bilateral filtering. IEEE Trans. Image Process. 21, 4361–4368 (2012)
  * Guo et al. [2013] Guo, R., Dai, Q., Hoiem, D.: Paired regions for shadow detection and removal. IEEE Trans. Pattern Anal. Mach. Intell. 35, 2956–2967 (2013)
  * Yu et al. [2017] Yu, X., Li, G., Ying, Z., Guo, X.: A new shadow removal method using color-lines. In: Computer Analysis of Images and Patterns, pp. 307–309 (2017)
  * Sasi and Govindan [2016] Sasi, R.K., Govindan, V.K.: Shadow Removal Using Sparse Representation Over Local Dictionaries. Engineering Science and Technology, an International Journal 19(2), 1067–1075 (2016)
  * Fan et al. [2019] Fan, H., Han, M., Li, J.: Image shadow removal using end-to-end deep convolutional neural networks. Applied Sciences (2019)
  * Qu et al. [2017] Qu, L., Tian, J., He, S., Tang, Y., Lau, R.W.H.: Deshadownet: A multi-context embedding deep network for shadow removal. In: Computer Vision and Pattern Recognition, pp. 2308–2316 (2017)
  * Wang et al. [2018] Wang, J., Li, X., Yang, J.: Stacked conditional generative adversarial networks for jointly learning shadow detection and shadow removal. In: Computer Vision and Pattern Recognition, pp. 1788–1797 (2018)
  * Isola et al. [2017] Isola, P., Zhu, J.-Y., Zhou, T., Efros, A.A.: Image-to-image translation with conditional adversarial networks. In: Computer Vision and Pattern Recognition, pp. 5967–5976 (2017)
  * Malherbe and Dalmasse [2019] Malherbe, J.-M., Dalmasse, K.: The new 2018 version of the meudon spectroheliograph. Solar Physics 294 (2019)
  * [15] pngkey.com: Transparent Cloud Texture Download. https://www.pngkey.com/maxpic/u2a9o0u2i1i1y3r5/
  * Ronneberger et al. [2015] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: Navab, N., Hornegger, J., Wells, W.M., Frangi, A.F. (eds.) Medical Image Computing and Computer-Assisted Intervention (MICCAI) (2015)
  * Heusel et al. [2017] Heusel, M., Ramsauer, H., Unterthiner, T., Nessler, B., Hochreiter, S.: Gans trained by a two time-scale update rule converge to a local nash equilibrium. Advances in neural information processing systems 30 (2017)


# WildGEN: Long-horizon Trajectory Generation for Wildlife

Ali Al-Lawati<EMAIL_ADDRESS>1234-5678-9012 Penn StateP.O. Box 1212State
CollegePennsylvaniaUSA16802 ,  Elsayed Eshra<EMAIL_ADDRESS>Penn StateP.O.
Box 1212State CollegePennsylvaniaUSA16802 and  Prasenjit Mitra<EMAIL_ADDRESS>Penn StateP.O. Box 1212State CollegePennsylvaniaUSA16802

(2023)

###### Abstract.

Trajectory generation is an important concern in pedestrian, vehicle, and
wildlife movement studies. Generated trajectories help enrich the training
corpus in relation to deep learning applications, and may be used to
facilitate simulation tasks. This is especially significant in the wildlife
domain, where the cost of obtaining additional real data can be prohibitively
expensive, time-consuming, and bear ethical considerations. In this paper, we
introduce WildGEN: a conceptual framework that addresses this challenge by
employing a Variational Auto-encoders (VAEs) based method for the acquisition
of movement characteristics exhibited by wild geese over a long horizon using
a sparse set of truth samples. A subsequent post-processing step of the
generated trajectories is performed based on smoothing filters to reduce
excessive wandering. Our evaluation is conducted through visual inspection and
the computation of the Hausdorff distance between the generated and real
trajectories. In addition, we utilize the Pearson Correlation Coefficient as a
way to measure how realistic the trajectories are based on the similarity of
clusters evaluated on the generated and real trajectories.

Trajectory Generation, Wildlife Trajectory, Trajectory Framework

††copyright: acmcopyright††journalyear: 2023††doi: XXXXXXX.XXXXXXX††price:
15.00††isbn: 978-1-4503-XXXX-X/18/06

##  1\. Introduction

Animals in the wild move for a variety of reasons, including foraging,
exploration, or migration (Venter et al., 2015). A better understanding of the
drivers of wildlife movement can have tremendous benefits in environmental
planning, preservation efforts, and anti-poaching. It contributes to
mitigating conflict between wildlife and human settlements (Main et al.,
2020).

Scientists in wildlife studies use datasets of animal trajectories collected,
typically, using telemetry devices, such as GPS collars, to capture the
animals’ exact location over time. However, the limitations of using these
devices arise from the complications and dangers associated with attaching the
devices to the animals, including constraints on device size and weight
depending on the animal, potential failures, and battery life.

As a result, considerable work in movement ecology attempts to generate
wildlife trajectories based on concepts such as Correlated Random Walk (CRW)
(Renshaw and Henderson, 1981) and Brownian Motion (Horne et al., 2007).
However, some of these methods are limited in their effectiveness to land-
based animals. Furthermore, in their basic form, these methods can only model
the wandering effect of animals within a foraging area (Bailey and Thompson,
2006). Methods such as Brownian Bridges (Technitis et al., 2015) extend these
movement models using a beginning and endpoint, and interpolate the in-between
points based on the maximum speed and time budget of the animal over the
expected duration of the trajectory. However, these methods fail to take
existing measured trajectories as input in the process of trajectory
generation. Instead, they rely on random variables to induce movement, and
thereby achieve limited success.

On the other hand, pedestrian and vehicle trajectory generation problems are
constrained by the road and pedestrian networks which are not applicable in
the wildlife domain. Likewise, trajectory generation in the UAV literature is
irrelevant as it is more concerned with optimization, and minimizing flying
time through the target waypoints (Singh and Fuller, 2001).

Figure 1. Real (green) vs Generated (blue) Trajectories

In this paper, we present a novel framework we call WildGEN that attempts to
generate trajectories by learning the statistical properties of a sparse set
of known trajectories. The framework, WildGEN, is based on proven deep
learning-based Variational Auto-encoders (VAEs) which have successfully been
applied to similar problems (Chen et al., 2021). WildGEN applies a post-
processing filter on each trajectory path using a Minimum Bounding Region
(MBR) based constraint. Furthermore, in order to reduce the excessive
wandering effect, a smoothing filter is applied.

In the next section, we discuss some related work, followed by a brief
motivation of the methods we used to develop our Framework. Section 4 presents
the experimental settings, experimental results, and an ablation study that
demonstrates the contribution of each component to the overall results. We
conclude and discuss future work in Section 5.

The contributions of this work can be summarized as follows:

  * •

We propose a framework that generates realistic long-horizon wildlife
trajectories based on a sparse set of real trajectories.

  * •

We experimentally test our framework and we validate the empirical results
using Hausdorff distance and the Pearson Correlation Coefficient metrics to
measure path and cluster similarities, respectively between the real and
generated trajectory sets.

  * •

We benchmark our method against Levy Walk/Flight (Zaburdaev et al., 2015),
which is state-of-the-art in long-horizon trajectory generation in the
wildlife domain. We also evaluate our methods over kriging (Iglesias et al.,
2018), using a Heteroscedastic Gaussian Process Regression (GPR) of multiple
trajectories on the same time horizon.

Figure 1 visually demonstrates the potential of this approach. We further
evaluate our approach using well-established trajectory similarity measures.

###  1.1. Problem Statement

Given a set of $n$ real trajectories, each consisting of $m$ $(x,y)$ points,
the objective is to produce a set of synthetic trajectories that maximize the
following:

  1. (1)

(Path Similarity) How similar are the points of the synthetic trajectories to
the real trajectories?

  2. (2)

(Cluster Similarity) How well do the synthetic trajectories align with the
real trajectories in terms of areas of concentration?

Formally, given $n$ real trajectories, denoted as

$T_{1},T_{2},\ldots,T_{n}$

where each $T_{i}$ consists of $m$ points:

$T_{i}=(p_{1}^{i},p_{2}^{i},\ldots,p_{m}^{i})\\}$

for $1\leq i\leq n$, and $p_{j}^{i}(j\in\\{1,2,\ldots,m\\})$ is a spatio-
temporal point consisting of latitude and longitude values. Find a set of:

$S_{1},S_{2},\ldots,S_{n}$

such that we maximize the similarity between the two sets:

$\max\text{ Similarity}(T,S)$

where $T=\\{T_{1},T_{2},\ldots,T_{n}\\},S=\\{S_{1},S_{2},\ldots,S_{n}\\}$, and
similarity encompasses Path and Cluster similarity. For the sake of
simplicity, we assume the number of synthetic trajectories is equal to the
number of real trajectories. We discuss the metrics used in Section 3.5.

##  2\. Related Works

Trajectory generation in the wildlife domain is dominated by models based on
Correlated Random Walk (CRW) (Bergman et al., 2000), Brownian motion (Horne et
al., 2007), and Levy Walk/Flight (Wosniack et al., 2017). These methods are
based on random steps and turning angles that follow a pre-defined
distribution of step length. Trajr is a popular wildlife trajectory R library
that provides tools to generate trajectories using each of the aforementioned
models (McLean and Skowron Volponi, 2018).

In (Technitis et al., 2015), the authors proposed a CRW trajectory generation
method based on the maximum speed and expected direction of movement. However,
the approach requires known starting and endpoints and is based on random
distributions of velocity and movement angle. Nonetheless, such methods are
more effective for simpler concerns such as missing data imputation, i.e. the
estimation of missing points between two known points. They are more useful in
capturing the micro movement of animals within their home range (Reynolds,
2009).

On the other hand, the generation of trajectory in the pedestrian literature
is constrained by road, other agents, and physical constraints. In (Chen et
al., 2021), an additional embedding layer is introduced to the VAE Encoder to
account for the road network. This ensures the road features are embedded in
the learning and understood by the generated trajectories. This has limited
applicability in wildlife movement where constraints are terrain-based. On the
other hand, a minimum bounding region (MBR) is capable of establishing
constraints on the spatial areas that a trajectory may be allowed to navigate.

Generative Adversarial Networks (GANs) are another class of generative
algorithms that have been utilized for the trajectory generation problem (Wang
et al., 2021). However, GANs are well known for their poor convergence
properties (Zhang et al., 2018). In addition, GANs generally require a large
sample of training data and are thus ineffective for a sparse sample set (Chen
et al., 2021).

Figure 2. The proposed WildGEN Framework

:
---

##  3\. Methodology

Figure 2 depicts the architecture of the proposed WildGEN to generate and
post-process a set of real trajectories. The framework is composed of four
main components:

  * •

A Variational Auto Encoder (VAE): takes a data set consisting of a training
set of n-day long trajectories, each of which represents a single
latitude/longitude pair per day for each subject.

  * •

A Gaussian Mixture Model (GMM): learns the probability distribution of the
training set in the latent space.

  * •

A trajectory smoothing filter: utilizes Savitzky-Golay smoothing to reduce the
excessive wandering effect observed on the generated trajectories.

  * •

A Minimum Bounding Region (MBR): analyzes the generated trajectories violation
of the constraints on acceptable movement region, and excludes outliers.

The framework consists of a training phase where the embeddings of each
trajectory are learned in the latent space. This allows the GMM to learn the
distribution of the training set in the latent space, and hence generate new
trajectories based on the distribution. The trajectories are decoded using the
learned weights from the decoding space. A smoothing post-processing step is
utilized on each synthetic trajectory, followed by an MBR step that accepts or
rejects it.

###  3.1. Variational Autoencoders

Variational Autoencoders (VAEs) (Kingma et al., 2019) are a class of
generative models that have gained significant traction in the field of
trajectory generation due to their ability to capture complex data
distributions and generate new samples. VAEs consist of two primary
components: the Encoder and the Decoder.

####  3.1.1. Encoder

The Encoder, denoted as $q(z|T)$, takes an input trajectory $T$ and maps it to
a latent space representation $z$. The Encoder learns to capture the essential
features and patterns present in the trajectory data. It compresses the high-
dimensional trajectory space into a lower-dimensional latent space,
effectively creating a data-driven embedding of the input trajectory.

####  3.1.2. Decoder

The Decoder, denoted as $p(T|z)$, takes a point $z$ from the latent space and
generates a trajectory $T$ in the original data space. The Decoder’s role is
to reconstruct trajectories from the latent space representations, effectively
learning the weights that convert latent data back to the original data space.

Real (green) vs Generated (blue)

(a) *

Real (green) vs Generated (blue)

(b) *

Real (green) vs Generated (blue)

(c) *

[.33] [.33] [.33]

Figure 3.  Three different sets of Generated Trajectories mapped against a
background of Real Trajectories

###  3.2. Savitzky-Golay Smoothing

Trajectory data, often obtained from GPS or motion tracking systems, can be
inherently noisy due to measurement errors and environmental factors.
Savitzky-Golay smoothing is a valuable tool for enhancing trajectory data
quality while preserving essential trajectory characteristics (Schafer, 2011).

By applying Savitzky-Golay smoothing to trajectory data, one can effectively
reduce the noise present in position ($x_{i}$) or velocity ($v_{i}$)
measurements. This noise reduction is achieved through a convolution operation
using a set of smoothing coefficients ($c_{j}$) within a local window:

$y_{i}=\sum_{j=-n}^{n}c_{j}x_{i+j}$

where $y_{i}$ represents the smoothed trajectory data at point $i$, $x_{i+j}$
are the neighboring data points within a window of size $2n+1$, and $c_{j}$
are the coefficients of the polynomial of order $m$ that minimize the mean
squared error between the smoothed data and the noisy trajectory data.

The choice of polynomial order ($m$) and the window size ($2n+1$) in the
context of trajectory data should be carefully considered. Too aggressive
smoothing (high $m$) may lead to the loss of important trajectory details,
while too little smoothing (low $m$) may not adequately mitigate noise.
Savitzky-Golay smoothing, when appropriately tuned, is a valuable mathematical
tool for improving the interpretability and reliability of trajectory data.

###  3.3. Minimum Bounding Region

Minimum Bounding Rectangles are the smallest axis-aligned rectangles (in two
dimensions) or hyperrectangles (in higher dimensions) that enclose a set of
trajectory points. In this paper, we generalize them using a convex hull and
refer to them as Minimum Bounding Regions (MBR). By encapsulating trajectories
within MBRs, we establish constraints on movement, such as maximum and minimum
bounds in both spatial dimensions.

The advantages of using MBRs lie in their simplicity and efficiency. They
provide a compact representation of the potential spatial footprint of
trajectory data, making it easier to detect outliers or unusual behavior.
Moreover, MBRs facilitate quick assessments of whether trajectories remain
within predefined geographical bounds, aiding in the enforcement of
constraints in the environment, whether it is regional bounds, rivers,
mountain ranges, and other physical constraints.

###  3.4. Baselines

In order to assess the performance of our framework we compare the results
obtained to existing methods for animal trajectory generation. As discussed
earlier, CRW-based and Brownian motion methods fail at generating long-horizon
trajectories for tasks such as migration. Similarly, deep learning methods
require large training samples.

As such we utilize two baselines that appear to be plausible for the task of
long-horizon trajectory generation, namely:

  * •

Levy Walk/Flight

  * •

Heteroscedastic Multi-variate Gaussian Process Regression (Heteroscedastic
GPR)

Next, we discuss the details of each of these methods, accordingly.

####  3.4.1. Levy Walk/Flight

is a stochastic process that models movement using a Cauchy distribution to
generate trajectory points, which are determined by a combination of step
length, variance, and angular standard deviation. These parameters are derived
from empirical data or real-world samples and are crucial inputs for the Levy
trajectory generation process.

The algorithm combines these parameter values of step length, variance, and
angular standard deviation with actual data points to produce the Levy
trajectory, which is a dynamic and realistic path that captures the
distinctive properties of the moving animal.

####  3.4.2. Heteroscedastic GPR

Heteroscedastic GPR (Le et al., 2005) is an extension of traditional Gaussian
Process Regression (GPR) that explicitly models and accounts for varying
levels of noise in the data. GPR is a powerful and versatile non-parametric
method commonly used for modeling complex relationships in data.
Traditionally, GPR assumes that the noise in the observed data is
homoscedastic, meaning that the variance of the noise is constant across all
data points. However, real-world data often exhibits heteroscedasticity, where
the variance of the noise varies across the input space or with respect to the
magnitude of the predictions. In such cases, traditional GPR may provide
suboptimal results, leading to biased predictions and inaccurate uncertainty
estimates.

Heteroscedastic GPR addresses the issue of varying noise levels by modeling
the noise variance as a function of input variables. This is achieved through
the introduction of a separate Gaussian process, often referred to as the
noise process or the noise model. The noise process provides an estimate of
the noise variance associated with each data point, which is then used to
adjust the likelihood function in the GPR model.

Given a dataset of $N$ observations $\\{(x_{i},y_{i})\\}_{i=1}^{N}$, where
$x_{i}$ represents the input variable and $y_{i}$ represents the observed
output variable. The observed data $y_{i}$ is affected by heteroscedastic
noise, meaning that the variance of the noise $\epsilon_{i}$ depends on the
input variable $x_{i}$.

The model can be characterized by:

  * •

Latent Function: The underlying latent function $f(x)$ is modeled as a
Gaussian process with a mean function $\mu(x)$ and a covariance (kernel)
function $k(x,x^{\prime})$, as in traditional GPR.

  * •

Noise Model: A separate Gaussian process is introduced to model the noise
variance as a function of the input variables. This is represented as $g(x)$,
where $g(x_{i})$ models the noise variance for each data point $x_{i}$.

The likelihood function in Heteroscedastic GPR is modified to account for the
varying noise levels:

$p(y_{i}|f(x_{i}),g(x_{i}))=\frac{1}{\sqrt{2\pi
g(x_{i})}}\exp\left(-\frac{(y_{i}-f(x_{i}))^{2}}{2g(x_{i})}\right)$

Here:

  * •

$p(y_{i}|f(x_{i}),g(x_{i}))$ is the likelihood of observing $y_{i}$ given the
latent function $f(x_{i})$ and the noise variance $g(x_{i})$.

  * •

$g(x_{i})$ represents the estimated noise variance for data point $x_{i}$.

  * •

$f(x_{i})$ is the value of the latent function at $x_{i}$.

  * •

$y_{i}$ is the observed data point at $x_{i}$.

The need for Heteroscedastic GPR arises from the requirement to model multiple
data points for a given input, where the input ($x_{i}$) is the day number in
the trajectory, and the output ($y_{i}$) is the set of all real points from
that day.

###  3.5. Evaluation Metrics

We base our selection of evaluation based on the objectives of trajectory
generation, which we reiterate here for convenience:

  1. (1)

(Path Similarity) How similar are the points of the synthetic trajectories to
the real trajectories?

  2. (2)

(Cluster Similarity) How well do the synthetic trajectories align with the
real trajectories in terms of areas of concentration?

In order to measure the first point, we use Hausdorff Distance similar to
(Chen et al., 2021), and we utilize the Pearson Correlation Coefficient to
measure the similarity of the cluster sets for the second metric. Next, we
describe both of them in further detail.

####  3.5.1. Hausdorff Distance

Hausdorff Distance is a valuable metric for assessing the similarity or
dissimilarity between movement paths. Trajectory data can be compared using
the Hausdorff distance to understand how closely two trajectories align or
diverge. The Hausdorff distance thus provides a way to measure the
dissimilarity between two trajectories by considering the maximum minimum
distance between points of one trajectory to the other trajectory.

The mathematical formulation of the Hausdorff distance (Hausdorff metric)
between two sets of points, A and B, in a metric space can be expressed as
follows:

Given two sets, A and B, each containing a number of points the Hausdorff
distance from set A to set B, denoted as H(A, B), is defined as:

Hausdorff distance from set $A$ to set $B$:

$H(A,B)=\max_{a\in A}\left(\min_{b\in B}d(a,b)\right)$

Hausdorff distance from set $B$ to set $A$:

$H(B,A)=\max_{b\in B}\left(\min_{a\in A}d(b,a)\right)$

where d(a,b) is the distance function that computes the distance between two
points.

####  3.5.2. Coefficient of Correlation

Pearson’s coefficient of correlation (Asuero et al., 2006), denoted as ”r” is
a fundamental statistical measure used to quantify the strength and direction
of linear relationships between two distributions. It is a widely used
statistical measure that quantifies the strength and direction of a linear
relationship between two distributions. Here are the key characteristics and
aspects of Pearson’s correlation coefficient:

  * •

An “$r$” value of -1 indicates a perfect negative linear relationship, where
one variable decreases as the other increases in a perfectly linear fashion.

  * •

An “$r$” value of 1 signifies a perfect positive linear relationship, where
both variables increase together linearly.

  * •

An “$r$” value of 0 implies no linear relationship, indicating that the
variables are independent of each other.

The calculation of Pearson’s correlation coefficient is based on the formula:

$r=\frac{\sum{(X_{i}-\bar{X})(Y_{i}-\bar{Y})}}{\sqrt{\sum{(X_{i}-\bar{X})^{2}}\sum{(Y_{i}-\bar{Y})^{2}}}}$

Here, $X_{i}$ and $Y_{i}$ represent individual data points from the two
variables, and $\bar{X}$ and $\bar{Y}$ are their respective means.

##  4\. Experiments

###  4.1. Datasets

The dataset used for our validation is The Western Palearctic greater white-
fronted geese dataset (Kölzsch et al., 2019). The dataset is hosted and was
obtained from Movebank111https://www.movebank.org. The dataset contains
reading of 91 different birds over a period of several months to multiple
years. We preprocessed the dataset using Moveapps (Kölzsch et al., 2022) to
limit the observations to one observation per day. Next, the time frames were
clipped to a window between March 1 to September 1 (185 days) to capture the
northward migration. Any subjects that had significant readings missing during
this time frame were dropped. Subjects that were tracked over multiple years
were split into multiple samples. Simple data interpolation was utilized to
fill in any remaining isolated gaps in the trajectories. The resultant dataset
consists of 60 trajectories over a period of 185 days.

###  4.2. Experimental Settings

Here we discuss the implementation details of the WildGEN framework as well as
the two baselines we compare it against.

Figure 4. Real (green) vs Generated (red) Trajectories using Levy Walk/Flight

####  4.2.1. Levy Walk/Flight Implementation Details

Figure 5. A Histogram of the step length distribution. Note that both axis
have a power multiple

A simple way of implementing Levy Walk/Flight involves analyzing the step
lengths of the trajectory and learning the coefficient of Cauchy distribution
($\alpha$) from that data. However, the $\alpha$ value computed for the
dataset was less than 1, which suggests that your step lengths have a very
heavy-tailed distribution with an increased likelihood of extremely long steps
as depicted by Figure 5. Nonetheless, for the purposes of this work, we used
the Trajr R library (McLean and Skowron Volponi, 2018), fed with the following
input variables:

  * •

Step Length

  * •

Angular Error

  * •

Linear Error

These variables were measured based on the training dataset.

As per the documentation of Trajr (McLean and Skowron Volponi, 2018), the
generated trajectory is computed with direction = 0. In order to guide the
movement in the correct northeast direction, the azimuth angle between the
average starting point and the average ending point was computed. The
generated trajectories were rotated counterclockwise with the computed angle.

Similar to WildGEN, the post-processing steps including Smoothing and MBR
inclusiveness check were applied to the generated trajectories, and
trajectories that violated the MBR were discarded.

Figure 4 demonstrates the generated trajectories from Levy Walk/Flight.

Figure 6. Real (green) vs Generated (red) Trajectories using Heteroscedastic
GPR

####  4.2.2. Heteroscedastic GPR Implementation Details

The Python GPy222https://github.com/SheffieldML/GPy library was used to model
Heteroscedastic GPR using the $GPHeteroscedasticRegression$ class with a
trained MLP and Bias kernels optimized on the training data. The GPR utilized
a random set of 50% of the training data. This reduced the computational
requirements needed to train the Heteroscedastic GPR.

In addition, the variances were computed using the norm of the $lng$ and $lat$
at each point. This helped limit overfitting of the model.

It is noteworthy that we also attempted to use the mean trajectory to model
the GPR and train a Matern kernel based on the mean squared error against each
training data sample. This approach was overfitting, and the samples generated
were almost indistinguishable.

Figure 6 depicts trajectories generated by Heteroscedastic GPR compared to
real trajectories. It can be observed that the clusters formed do not visually
appear to align with the clusters of the real trajectories.

Figure 7. Training MSE Figure 8. Compressed Space

####  4.2.3. WildGEN Implementation Details

In the context of this study, the input structure comprises a dataset
consisting of 60 records, each encompassing longitude and latitude
coordinates, spanning a period of 185 days. These spatial-temporal records are
meticulously rearranged into a vector representation composed of 370 inputs.
It is noteworthy that prior to further processing, the input data undergoes a
normalization step, which involves dividing each input value by 0.3 times the
maximum input value, ensuring that the data adheres to a standardized scale.

The neural network architecture takes the form of an VAE. It initiates with an
input layer of 370 neurons, which serves as the point of entry for the data.
Sequentially, the network progresses through a hidden layer featuring 300
neurons, followed by another hidden layer consisting of 100 neurons.
Subsequently, the architecture transitions into a lower-dimensional latent
space representation, characterized by a dimensionality of 3. Following the
latent space encoding, the network proceeds with two additional hidden layers,
each comprising 100 neurons, facilitating the extraction of essential data
features. The network culminates in a final hidden layer containing 300
neurons, serving as the penultimate stage. Ultimately, the network’s journey
concludes with an output layer designed to match the initial input’s
dimensionality, comprising 370 units. This VAE architecture effectively
captures data patterns, allowing for both data compression and reconstruction.

Activation functions are carefully chosen for each layer to regulate the flow
of information through the neural network. The first hidden layer is governed
by a modified Rectified Linear Unit (ReLU) activation function, while the
second hidden layer employs a linear activation function. The latent space is
also characterized by a linear activation function. Notably, the fourth and
fifth hidden layers utilize a modified ReLU activation function, which is
characterized in this work by a slope of 0.06 for positive values of the
argument and a slope of 0.001 for negative values.

In this VAE-based study, the Mean Squared Error (MSE) is the primary loss
function, assessing data reconstruction quality. The MSE, as shown in Figure7,
quantifies the model’s ability to successfully replicate input data within its
latent space. Additionally, a Gaussian Mixture Model (GMM) with 15 components
is fitted to the latent space using the 60 compressed samples shown in Figure
8, for deeper insight into data distribution. To generate trajectories, we
employ a sampling approach within the latent space, drawing samples from the
GMM, as showcased in Figure 8, and subsequently propagate these samples
through the decoder network.

Smoothing: The Savitzky-Golay Smoothing was implemented against each
trajectory based on a manual visual tuning of the parameters, namely window
size and polynomial order. As Fig 9 demonstrates, the filter can reduce
excessive wandering of the generated trajectories in comparison with the real
trajectories.

Figure 9. Savitzky-Golay Smoothing Filter: Before (blue) and After (red) for a
sample generated trajectory

MBR: The Minimum Bounding Region (MBR) was calculated from the real
trajectories to constrain the region in which the birds can fly. In different
problems, one or multiple MBRs can be used to signify water bodies,
impenetrable mountain ranges, and other terrain-based obstacles. The MBR step
discards any trajectory that does not fully lie inside the calculated MBR.
This has resulted in discarding 29.6% of the generated trajectories and
retaining the remaining trajectories.

Figure 10. A Minimum Bounding Region based on the real trajectories

:
---

###  4.3. Model Comparison

(a)

(b)

Figure 11. Choosing the best value of K for K-means clustering based on two
separate methods

In order to compare the models, the first 60 eligible trajectories from each
model’s synthetic set were used to match the set of real trajectories on the
evaluation metrics. The remaining trajectories were not used during the
benchmarking process in this step.

####  4.3.1. Path Similarity

The Hausdorff distance between each pair of trajectories from the synthetic
and real set were compared. Only the lowest value was retained for each
synthetic trajectory. In other words, for each generated trajectory, the
distance was computed with the trajectory closest to it from the real
trajectory set.

We report the minimum, mean, and average distance for each model in Table 2.

####  4.3.2. Cluster Similarity Calculation

K-means clustering was performed on the set of all points, i.e.

$\bigcup_{i}\\{p_{1}^{i},p_{2}^{i},\ldots,p_{n}^{i}\\}$

for all $T_{i}$. The value of $K$ was chosen using a Silhouette Score and
validated visually using the Elbow method (Saputra et al., 2020). The number
of clusters based on the geographic projection of the real trajectory should
be between 10 and 15; K=13 was used in our experiments.

The trained clustering model was applied to the generated trajectory points of
each model, and the Pearson Correlation Coefficient was calculated based on
the distribution of points in each cluster.

####  4.3.3. Results

From Table 2, we can observe that WildGEN provides a better overall Hausdorff
Distance on average than the other models. Further, the clusters of WildGEN
have better correspondence to the real trajectories. While Heteroscedastic GPR
performs well on the Hausdorff distance metric, the results are influenced by
the low variance in the model. It is visually apparent that trajectories
generated by Heteroscedastic GPR follow the same path with little variance
through the first half of the trajectory as depicted in Figure 6.

On the other hand, it can be observed that Levy Walk/Flight fails to generate
any realistic trajectories. The metrics also indicate that the trajectories
generated fail to capture any properties of the real trajectories. This is not
surprising given the limited amount of information about the real trajectories
that can be fed into the algorithm.

However, WildGEN performs well on the metrics in line with visual inspection.
It performs significantly better than Levy Walk/Flight on average Hausdorff
Distance and provides a 15.5% improvement on average compared to
Heteroscedastic GPR. There is a very high correlation between the cluster
distributions and the generated points are closer on average to real
trajectory points which is almost double the cluster similarity observed using
Heteroscedastic GPR.

Table 1. Test Results for Hausdorff Distance and Pearson Correlation Coefficient Method | Hausdorff Distance | Pearson (r)
---|---|---
Min | Max | Avg |
Levy Walk/Flight | 22.251 | 39.799 | 30.590 | -0.0669
Heteroscedastic GPR | 5.760 | 7.155 | 6.208 | 0.4563
WildGEN | 3.542 | 7.450 | 5.244 | 0.9936

(a) Raw

(a) *

(b) Smoothed

(b) *

(c) MBR Filtered

(c) *

(d) WildGEN

(d) *

[.24] [.24] [.24] [.24]

Figure 12. The effect of different components on WildGEN on the generated
trajectories

###  4.4. Ablation Study

In the ablation study, we want to measure the effects of Smoothing and Minimum
Bounding Region (MBR) on the metrics measured.

Table 2. Test Results for Hausdorff Distance and Pearson Correlation Coefficient Method | Hausdorff Distance | Pearson (r)
---|---|---
Min | Max | Avg |
No Post Processing | 0.306 | 9.440 | 4.172 | 0.9881
Smoothing Only | 3.542 | 10.047 | 5.565 | 0.9865
MBR Only | 0.395 | 6.425 | 3.849 | 0.9928
WildGEN | 3.542 | 7.450 | 5.244 | 0.9936

In order to perform the ablation study correctly, we ran the experiments on a
random sample of trajectories, applied the specified method on the same set,
and retained only the first 60 trajectories in line with the experiments in
the previous section.

The outcomes of our ablation study reveal that smoothing alone yields
diminished results of the Hausdorff distance and Pearson Correlation
Coefficient. However, as discussed previously, smoothing effectively addresses
a specific visual challenge, notably the mitigation of excessive wandering. To
provide a mathematical evaluation of the advantages provided by smoothing, we
may need to consider additional metrics.

On introducing MBR to our trajectories, the measurable impact of smoothing
diminishes. Nevertheless, the combined application of both smoothing and MBR
yields a small improvement in the Pearson Correlation Coefficient. It can be
observed that the penalization against the Hausdorff distance metric extends
to the WildGEN scenario.

Our findings also serve to validate the viability of the raw trajectories
generated by the VAE, as corroborated by the metrics introduced. Notably, the
utilization of smoothing significantly contributes to their improved visual
inspection when compared to real trajectory data.

It is important to note that comparing Smoothing Only and MBR Only scenarios
is not completely objective. This is because MBR excludes many trajectories
that the smoothing step retains. On the other hand, if we doctor the random
set to ensure that only trajectories that meet the MBR criteria are used in
the ablation study, it would defeat the purpose of the improvements of MBR. In
that case, MBR and No Post Processing would have identical results.

##  5\. Conclusions and Future Work

This work presented a novel method by which synthetic trajectories may be
generated on a long horizon. The framework, WildGEN, can be characterized as
follows:

  * •

Generate near-real, long-horizon trajectories based on a relatively small
number of real trajectories

  * •

Performs post-processing steps to further improve the generated trajectories

  * •

Shows improved accuracy compared to the state-of-the-art methods

Based on the findings from this work, there are multiple avenues for future
work. A key limitation identified is the evaluation metrics used are not
sufficiently representative of a good trajectory. In particular, Hausdorff
distance fails to capture the enhancement in the trajectory, and it is
necessary to test other metrics such as Frechet Distance to better assess the
effect of smoothing. It is also interesting to identify and utilize other
potential smoothing filters.

In addition, as a follow-up to this work, we intend to:

  * •

Utilize other generative models, including normalizing flows that can provide
the pertinent statistical weights of the generated trajectories

  * •

Extend and further investigate the applicability of WildGEN for various
wildlife trajectory scenarios, including land-based and water-based animals.

## References

  * (1)
  * Asuero et al. (2006) Agustin Garcia Asuero, Ana Sayago, and AG González. 2006\.  The correlation coefficient: An overview.  _Critical reviews in analytical chemistry_ 36, 1 (2006), 41–59.
  * Bailey and Thompson (2006) Helen Bailey and Paul Thompson. 2006.  Quantitative analysis of bottlenose dolphin movement patterns and their relationship with foraging.  _Journal of Animal Ecology_ 75, 2 (2006), 456–465.
  * Bergman et al. (2000) Carita M Bergman, James A Schaefer, and SN Luttich. 2000\.  Caribou movement as a correlated random walk.  _Oecologia_ 123, 3 (2000), 364–374.
  * Chen et al. (2021) Xinyu Chen, Jiajie Xu, Rui Zhou, Wei Chen, Junhua Fang, and Chengfei Liu. 2021\.  TrajVAE: A Variational AutoEncoder model for trajectory generation.  _Neurocomputing_ 428 (2021), 332–339.  https://doi.org/10.1016/j.neucom.2020.03.120
  * Horne et al. (2007) Jon S Horne, Edward O Garton, Stephen M Krone, and Jesse S Lewis. 2007.  Analyzing animal movements using Brownian bridges.  _Ecology_ 88, 9 (2007), 2354–2363.
  * Iglesias et al. (2018) I Iglesias, F Montes, M Martínez, A Perez, A Gogin, D Kolbasov, and A de la Torre. 2018.  Spatio-temporal kriging analysis to identify the role of wild boar in the spread of African swine fever in the Russian Federation.  _Spatial statistics_ 28 (2018), 226–235.
  * Kingma et al. (2019) Diederik P Kingma, Max Welling, et al. 2019\.  An introduction to variational autoencoders.  _Foundations and Trends® in Machine Learning_ 12, 4 (2019), 307–392.
  * Kölzsch et al. (2022) Andrea Kölzsch, Sarah C Davidson, Dominik Gauggel, Clemens Hahn, Julian Hirt, Roland Kays, Ilona Lang, Ashley Lohr, Benedict Russell, Anne K Scharf, et al. 2022\.  MoveApps: a serverless no-code analysis platform for animal tracking data.  _Movement Ecology_ 10, 1 (2022), 30.
  * Kölzsch et al. (2019) A Kölzsch, GJDM Müskens, S Moonen, H Kruckenberg, P Glazov, and M Wikelski. 2019\.  Data from: Flyway connectivity and exchange primarily driven by moult migration in geese [North Sea population].  https://doi.org/doi:10.5441/001/1.ct72m82n
  * Le et al. (2005) Quoc V Le, Alex J Smola, and Stéphane Canu. 2005.  Heteroscedastic Gaussian process regression. In _Proceedings of the 22nd international conference on Machine learning_. 489–496.
  * Main et al. (2020) Michael T. Main, Robert A. Davis, David Blake, Harriet Mills, and Tim S. Doherty. 2020.  Human impact overrides bioclimatic drivers of red fox home range size globally.  _Diversity and Distributions_ 26, 9 (2020), 1083–1092.  Publisher: Wiley Online Library.
  * McLean and Skowron Volponi (2018) Donald James McLean and Marta A Skowron Volponi. 2018.  trajr: an R package for characterisation of animal trajectories.  _Ethology_ 124, 6 (2018), 440–448.
  * Renshaw and Henderson (1981) Eric Renshaw and Robin Henderson. 1981.  The correlated random walk.  _Journal of Applied Probability_ 18, 2 (1981), 403–414.
  * Reynolds (2009) A. M. Reynolds. 2009\.  Scale-free animal movement patterns: Lévy walks outperform fractional Brownian motions and fractional Lévy motions in random search scenarios.  _Journal of Physics A: Mathematical and Theoretical_ 42, 43 (Oct. 2009), 434006.  https://doi.org/10.1088/1751-8113/42/43/434006
  * Saputra et al. (2020) Danny Matthew Saputra, Daniel Saputra, and Liniyanti D Oswari. 2020.  Effect of distance metrics in determining k-value in k-means clustering using elbow and silhouette method. In _Sriwijaya International Conference on Information Technology and Its Applications (SICONIAN 2019)_. Atlantis Press, 341–346.
  * Schafer (2011) Ronald W Schafer. 2011\.  What is a Savitzky-Golay filter?[lecture notes].  _IEEE Signal processing magazine_ 28, 4 (2011), 111–117.
  * Singh and Fuller (2001) Leena Singh and James Fuller. 2001.  Trajectory generation for a UAV in urban terrain, using nonlinear MPC. In _Proceedings of the 2001 American Control Conference.(Cat. No. 01CH37148)_ , Vol. 3. IEEE, 2301–2308.
  * Technitis et al. (2015) Georgios Technitis, Walied Othman, Kamran Safi, and Robert Weibel. 2015.  From A to B, randomly: a point-to-point random trajectory generator for animal movement.  _International Journal of Geographical Information Science_ 29, 6 (2015), 912–934.
  * Venter et al. (2015) Jan A. Venter, Herbert HT Prins, Alla Mashanova, Willem F. de Boer, and Rob Slotow. 2015.  Intrinsic and extrinsic factors influencing large African herbivore movements.  _Ecological Informatics_ 30 (2015), 257–262.  Publisher: Elsevier.
  * Wang et al. (2021) Xingrui Wang, Xinyu Liu, Ziteng Lu, and Hanfang Yang. 2021\.  Large scale GPS trajectory generation using map based on two stage GAN.  _Journal of Data Science_ 19, 1 (2021), 126–141.
  * Wosniack et al. (2017) Marina E Wosniack, Marcos C Santos, Ernesto P Raposo, Gandhi M Viswanathan, and Marcos GE Da Luz. 2017.  The evolutionary origins of Lévy walk foraging.  _PLoS computational biology_ 13, 10 (2017), e1005774.
  * Zaburdaev et al. (2015) Vasily Zaburdaev, Sergey Denisov, and Joseph Klafter. 2015\.  Lévy walks.  _Reviews of Modern Physics_ 87, 2 (2015), 483\.
  * Zhang et al. (2018) Zhaoyu Zhang, Mengyan Li, and Jun Yu. 2018.  On the convergence and mode collapse of GAN.  In _SIGGRAPH Asia 2018 Technical Briefs_. 1–4.