# Deep Learning based Automated Forest Health Diagnosis from Aerial Images

Chia-Yen Chiang1, Chloe Barnes2, Plamen Angelov1, and Richard Jiang1

(1LIRA, Lancaster University, Lancaster, LA1 4WY, UK
22Excel geo, Hall Farm 2, Sywell Aerodrome, Sywell, NN6 0BN)

###### Abstract

Global climate change has had a drastic impact on our environment. Previous
study showed that pest disaster occured from global climate change may cause a
tremendous number of trees died and they inevitably became a factor of forest
fire. An important portent of the forest fire is the condition of forests.
Aerial image-based forest analysis can give an early detection of dead trees
and living trees. In this paper, we applied a synthetic method to enlarge
imagery dataset and present a new framework for automated dead tree detection
from aerial images using a re-trained Mask RCNN (Mask Region-based
Convolutional Neural Network) approach, with a transfer learning scheme. We
apply our framework to our aerial imagery datasets,and compare eight fine-
tuned models. The mean average precision score (mAP) for the best of these
models reaches 54%. Following the automated detection, we are able to
automatically produce and calculate number of dead tree masks to label the
dead trees in an image, as an indicator of forest health that could be linked
to the causal analysis of environmental changes and the predictive likelihood
of forest fire.

Keywords: Deep learning, object detection, object segmentation, forestry,
biodiversity

##  1 Introduction

More and more evidence indicates that global climate change has led to
deforestation and extinction of species over the past decade. The incremental
changes in temperature have not only caused extreme events such as floods,
drought and gales, but inevitably affect landscapes and whole ecosystems. In
such a situation, species may struggle to adapt to changes [1]. Moreover, it
is reported [2] that 110,000 hectares of forest have died in 2018, which is
around the size as 150,000 football pitches or one nineteenth of Wales, due to
extreme heat and storms.

Certain species whose habitats were living woods have lost their habitats due
to deforestation and dying trees. On the other hand, organism or creatures
that take dead wood as habitats (for example: fungus, bacteria, predatory
beetles, and parasitic wasps) increase their population. Pests such as bark
beetles and cardinal beetle can nest in dead tree trunks, and then may move to
healthy trees and kill them [2].

In addition, the dry, inflammable nature of dead wood makes fire disasters
happen more frequently [3]. In California, poor management of dead wood in
forests has been identified as the major cause of wildfire. Both brought huge
costs since the fire destroyed neighboring healthy trees and spread fast to
villages nearby. Forest fires also trigger severe smog; for example, in 2015,
wildfire smoke in Indonesia caused an estimated 100,000 premature deaths [4].

Because the of impacts on ecosystems, inhabitants and households mentioned
above, technology to monitor forest areas has been developed rapidly in the
past decades. Technologies such as remote airborne sensor or satellite imagery
made forest monitoring less costly, and allows monitoring without human
involvement. By using these technologies, forest experts and disaster
specialists not only can monitor forests in real-time, but also adopt
preventive measures for natural disasters.

On the other hand, pest disaster, one of the natural disasters caused by
climate change, has recently drawn experts’ attention. The Canopy Health
Monitoring (CanHeMon) [5] project aimed to eliminate the danger to forests
from pine wood nematode (PWN), and so previous research has been done to
detect dead trees within forests, to prevent them from becoming vectors for
PWN. A MaxEnt-based iterative image analysis algorithm, applied to remote
sensing data, as able to detect individual declining tree crowns with high
accuracy [5]. The algorithm searches for the distribution with the maximum
entropy for a given input set of occurrences and a set of constraining
variables. As higher entropy indicates high disorder and abundance of
information, the algorithm can filter out irrelevant variables, leaving the
most useful variables highly correlated to dead trees. However, the results
revealed that the rate of false positives increases as the sample size
decreases. Moreover, eliminating small size samples (smaller than 2.5m2)
harmed the true positive detection rate, dropping from 80% to 65%. This may
indicate that the algorithm is not suitable to detect small objects or that
there were not enough samples of small trees in the training data.

In addition, the CanHeMon project suggested that using deep learning
algorithms may provide a more generalised image analysis model and allow us to
identify pixel-wise shapes of trees in an image. To investigate this, the
CanHeMon team also performed tests using a neural network model. The results
showed that rate of false positives was lower than for the MaxEnt-based
algorithm, but the behaviour of the loss function indicated more epochs and
iterations were needed in training the model. Hence, in this research, we
extend the deep learning methods applied by the previous team.

In this study, we used a Mask Region with Convolutional Neural Networks (Mask-
RCNN) algorithm, accompanied by transfer learning [6] to deal with the problem
of pixel-wise object identification and to enhance the learning pattern, i.e.
through additional training sessions, the overall score performance can be
improved by using pre-trained weights produced in the previous session.

Moreover, to cope with small datasets, and problems with unlabelled data, we
use synthetic tools: randomly assigning instances into backgrounds, both
extracted from raw data. At the same time, we output COCO format annotations
[7] to be used in inference. In general, manually annotating these masks is
quite time-consuming, frustrating and commercially expensive [8]. However, we
were inspired by Adam Kelly, the founder of Immersive Limit, and were able to
automatically create thousands of synthetic images and their annotations with
minimal manual effort. [9]
The study makes four major contributions:

(1) Use synthetic data to do Mask RCNN modeling and reach about 54% on mean
average precision score (mAP) . The results show that it is possible to use
synthetic methods to mitigate data scarcity.

(2) Mask prediction is limited by selected ROIs and can be found from high
loss of RPN. Therefore, applying semantic segmentation may be a more practical
way than instance segmentation to detect dead tree region.subsection 6.2

(3) Using two optimizers in training process caused loss converged faster than
using one optimizer. subsection 6.2

In addition, Figure 1 provides our framework for forest analysis. More details
will be in the following sections.

Figure 1: Structure of the Mask-RCNN model.

##  2 Related work

###  2.1 Satellite imaging for forest analysis

Although we did not use satellite images in our research, they have been used
to monitor forest cover for a long time. Analysing the change of forest cover
generally relies on visual interpretation of a dataset. Accompanied with groud
measurements, scientists are able to calculate the loss of forest cover over
time and doing assess properties from the satellite imagery, such as
biophysical and biochemical properties of forest. [10]

In addition, in various studies[11][12][13], Landsat imagery[14], a popular
open database of satellite imagery, was sampled or trained to monitor forest
cover changes. Landsat imagery resolution is between 15-60m with multiple
bands. ”Landsat is a joint effort of the U.S. Geological Survey (USGS) and the
National Aeronautics and Space Administration (NASA)…The data are useful to a
number of applications including forestry, agriculture, geology, regional
planning, and education. ”[15]

###  2.2 AI techniques in object detection

With the remarkable advancement in artificial intelligence techniques in the
past decades, there are several object detection algorithms which are now
generally used, for example R-CNN[16], Fast-RCNN[17], Faster-RCNN[18],
YOLO[19], SSD[20], and NAS-Net[21]. When performing object detection, we are
interested in finding two things: where our target objects are located, and
which categories they belong to[22]. Our goal is to draw bounding boxes around
target objects in an image, and output a percentage indicating how certain we
are that they belong to certain classes. Moreover, state of the art techniques
allow real-time prediction of object localisation and classification, for
example You Only Look Once (YOLO) and Single-Shot Detector (SSD).

However, the cost of the fast inference time of such methods may be lower
accuracy and low-resolution feature maps, due to a single network evaluation
[23].. They use a single stage to complete both object localisation and
classification, but in constrast, region based algorithms (RCNN, Fast-RCNN,
Faster-RCNN and Mask-RCNN) separate the process into to two stages: first
using a region proposal network (RPN) to propose a bounding box for a
candidate object; then predicting the class and box offset, and outputting a
binary mask for each region of interest’s ROI [24]. That means we calculate
both the classification scores and bounding box regression in second stage to
decide which ROIs have the highest probability to contain target objects, and
simultaneously perform instance segmentation on selected ROIs. Here, instance
segmentation [24] means to differentiate each target object and cut them apart
even when they are labeled in the same class. It is often compared to semantic
segmentation, labeling on each pixel with the class of its enclosing object or
region (see Figure 2). Note that we only took segmentation into account
because it was more precise to cut off dead tree regions in pixel-wise manner
than simply output bounding boxes; for example, some boxes may contain a big
proportion of living trees.

###  2.3 CNN and FCN in pixel-wise dead forest mapping

Recent studies [25][26][27][28]showed CNN and FCN were used in sementic
segmentation on forest mapping. In [25] a VGG neural network (an architecture
of CNN) was proposed to classify pixels with a sliding window technique. One
notable approach called ensemble learning was introduced to aggregate
different predictions generated by the same model. This produced determistic
estimates and reduced uncertainties caused by weight initialization settings.
The resulting global accuracy predicting health status of two tree species was
around 80% to 85%. However, dead trees were misclassified (defined as
”omission error”) at a rate of 35% to 40%.

Another study suggested sliding windows [29] cause high computational cost
because neighbour pixel information, which contains the same information as
centre pixels, is used to classify centre/target pixel. An algorithm such as
FCN is developed to prevent duplicate computation by predicting many pixels at
once. In addition, FCN can take various input sizes and output a dense
prediction.

A DenseNet (an architecture of FCN) proposed in [26] has outstanding results
for classifying dead trees and standing dead trees. The resulting recalls and
precisions are mainly within 90 to 100% although a recall of a subset was
67.8%. The author also encouraged that instance level prediction can be done
to improve future research.

###  2.4 Transfer learning

In the field of image processing, there have been many studies using transfer
learning. The technique uses an existing trained network to improve the
training accuracy for a new model.

Since an outstanding trained network such as COCO or ImageNet has been trained
using many images with a big number of labels, it seems sensible to reuse
fundamental parts of the pre-trained net as low-level feature extractors, i.e.
to recognize edges, colours and textures - whether detecting complex and
subtle objects such as a human fact or simple objects such as rectangles,
early layers in neural networks always detect the same basic shapes and edges.

To be more specific, the COCO dataset[7] is a large-scale object detection,
segmentation, and captioning dataset. It has over 200 000 labeled images with
over 80 labels, while the ImageNet dataset has tens of millions of images with
thousands of labels. [30]

[6] So to train a model on a custom dataset, all we need to do is replace the
classification and bounding box predictors with customised ones which are able
to predict our labels. First, we take a pre-trained net without its FC part,
since Fully Connected layers are used only for the old model’s classification.
Next, a new FC containing customized label classifiers replaces the old. Then,
a convolution base which is initially close to the old classifier needs to be
frozen since it contains high level features (for example cat ears, dog eyes,
or a more representative feature) related to classification in the old model
which may not be useful for our own customized classification problem.

Figure 2: Semantic segmentation (left) versus instance segmentation (right).
Image from [24].

###  2.5 Transfer learning in forest aerial imagery

We used a COCO or Imagenet pre-trained net and retrained the ”head” parts for
our model in the first session, i.e. we kept the COCO or Imagenet backbone but
retrained the RPN, classifier and mask heads of the network. [31]

In the first training session, we only trained the head since we want to
leverage the advantage from COCO and ImageNet, that is, use the weights that
have been trained from many images. However, we found that merely training the
head was not sufficient to get higher accuracy. Because of that, in the second
and third sessions, we train all layers twice. Therefore, we denote the deeper
model as “head + all + all”.

Then in the second training session, we used the last weight produced in the
first session to continue training. This gave us faster training times, and
allowed us to use a well-trained model even when we had scarce data. Related
studies on detecting objects in satellite imagery of urban landscapes[32] show
that by using an ImageNet pre-trained net, accuracy was increased by around 5%
for a fixed length of training time. However, the study also showed a less
satisfying performance even after using transfer learning because the model is
trained to predict objects from another dataset which has different domain
features, which may not be similar to our own training data. Therefore,
problems such as how transferable features in deep neural networks are, and
how to effectively reuse pre-trained nets when most features of the target
objects are different from those of the training data, have been discussed in
[33]. A Selective Learning Algorithm (SLA) has been proposed in this case,
which learns useful features from unlabelled data and introduces an
intermediate domain to close the gap between target objects and training data.

##  3 Preliminary on CNN and Mask RCNN

###  3.1 CNN

####  3.1.1 What is a convolutional neural network?

A convolutional neuron network (covnet, or CNN) is a class of deep neural
network that has proved very powerful for image classification. In 2012, Alex
Krizhevsky used CNN to win the ImageNet competition which made neural nets
come to prominence in the field of computer vision. [30]

####  3.1.2 Why use a convolutional operation?

We can view an image as having 3D inputs with (height, weight, channel) so the
input features have dimension H$\times$W$\times$C. For deeper layers, the
channel value is usually increased in order to learn more features. A channel
value is also the number of filters which contain special features we want to
detect. In addition, given pre-defined filters, we can obtain corresponding
feature maps by applying convolutional operation on input images. The
operation not only shrinks the input dimension but preserve the spatial
relationship between pixels; therefore a series of features are as well
preserved. Figure 3 shows different feature maps produced by applying filters
to the same raw images. Usually, a typical Convnet looks like Figure 4.

Figure 3: Different feature maps produced by applying various filters to the
same raw image. From [34]. Figure 4: Structure of a convolutional neural
network. From [35].

####  3.1.3 Why can’t we use a simple CNN to predict object location?

  * •

The number of target objects and their coordinates in the raw picture is
unknown before detection, so output nodes can’t be specified. In CNN, however,
a fixed number of nodes is predefined for the output layer. Normally, outputs
in CNN are fixed, predefined classes. Therefore, additional layers are needed
to locate their coordinates.

  * •

To classify if an object is a target or not, a number of fixed-sized input
images are fed to CNN. However, in the object detection problem, input objects
may have different sizes. In this case, a ROI pooling layer in RCNN is added
to scale all objects in fixed size before passing them to the classification
process.

###  3.2 Mask RCNN structure

In this section, we will first introduce Faster RCNN, the predecessor to Mask
RCNN; then we will discuss Mask RCNN. This will give a broad picture of how
Mask RCNN developed and the structure of the algorithm (see Figure 5).

####  3.2.1 Faster RCNN

Faster RCNN begins with a standard convolution neural network (CNN) for
feature extraction, which produces feature maps. A feature map is the output
activation for a given filter in each layer. When a filter matches certain
features on an image, the convolution operation [36] returns high activation
values and reveal those features in the feature maps. A feature map is also a
mapping of where a certain kind of feature is found in the image [34].

Then, it uses a Region Proposal Network (RPN) to extract Regions of Interest
(ROIs) and resize the ROIs to a fixed size using ROI pooling. (Note: In
section subsection 4.2 , we will talk about how RPN generates ROIs.)

Figure 5: Structure of a Mask RCNN model. From [24]. Figure 6: How a max
pooling layer works

A pooling layer then reduces the dimension of the ROIs and the number of
parameters in order to save training time and combat overfitting [35]. It
divides an ROI into a number of sections by pre-defined output size (Note: it
is not necessary that ROI is perfectly divisible by the section number.) Next,
a maximum value is found for each section and we replace all other pixel
values in the section with this value. [18] See Figure 6.

Finally, the ROIs after max-pooling are flatten and fed to two fully connected
layers (FCs). The flatten process makes pooled ROIs turn into vectors which
store non-spatial features for later calculate the multilinear relationship
between input neurons (in each ROI) and classes, and between neurons and
bounding box size.

One FC layer calculates a bounding box regressor to refine the bounding box
size, i.e. make each box covers the regions over its target object.
Coordinates information will be sent to bounding box loss caculation. In next
iteration, to minimize lost, bounding box coordinates will be updated
therefore produce more accurate bounding box; and the other FC layer predicts
a vector of probabilites over output classes using a softmax function. The
probabilies will be used in calculation of cross-entropy loss which we adjusts
probability information to minimize the cost in every iteration.

####  3.2.2 Mask-RCNN

Mask-RCNN can be thought of as a combination of a Faster RCNN that does object
detection (class + bounding box) and a Fully Convolutional Network (FCN), a
semantic segmentation method producing a pixel-wise mask. (See Figure 7)

In addition, the author of Mask RCNN [24] did some fine-tuning of Faster RCNN
and FCN. First, since the ROI pooling layer in Faster RCNN can lead to the
location of target objects in the feature map being misaligned from their
location in the original image, this is dealt with by refining the mask with
an ROI Align layer. ROI Align is a quantization-free layer, preserving exact
spatial locations by replacing quantization with bilinear interpolation [24].
Bilinear interpolation is used for computing the floating-point location
values in the input [33].

Then FCN, a semantic segmentation method, is applied to each ROI in order to
predict a pixel-wise, binary mask (see Figure 8). Unlike the original setting
of FCN, here the activation function in FCN was changed from softmax to
sigmoid which was shown to significantly increase Average Precision (AP) [24].
The reason is that FCN only calculates per-pixel sigmoid values and an average
binary cross-entropy loss for a single, wholly classified object each time,
rather than calculating per-pixel softmax and a multinomial loss for all
pixels. Because of this, label prediction is decoupled with mask generation
which avoids competition among classes.

Figure 7: Structure of the Mask-RCNN model. Image from [link]. Figure 8: FCN
predicts a pixel-wise binary mask. Image from [link].

####  3.2.3 Backbone of Mask-RCNN

Mask-RCNN uses Resnet-FPN as its default backbone and is a type of CNN which
produces feature maps; the structure of Feature pyramid networks (FPN) allows
us to build high-level semantic feature maps at different scales. When images
are fed into FPN, first they go into the bottom-up process: the process
functions the same as classic covnet extracts features, four feature maps (C2,
C3, C4, C5) are produced and get reduced size (1/2 the size of the previous
map) as the level goes up. At the same time, more features are extracted on
the upper level but have lower resolution because of the convolution operation
process.

In addition to the loss of resolution introduced by this, high resolution
feature maps from lower level are not used for detection because of the lack
of sufficient information (for example, it is possible that only some edges
are detected). These two disadvantages make detection of small objects very
difficult.

In the FPN paper [37], however, the author added a top-down pathway (see
Figure 9) to reconstruct a high resolution from a semantic rich layer by an
upsampling method. The method doubles the map size (twice the size of the
previous map) as the level goes down, with the resulting feature maps denoted
as (P2, P3, P4, P5). A bottom-up feature map which goes through a lateral
connection is then merged with a top-down feature map accordingly. The
connection simply applies a 1$\times$1 convolution filter on each bottom-up
feature map; therefore, it retains map sizes but reduces channel depths. In
addition, to deal with spatial aliasing (the resulting location of objects may
be different from their locations in the raw image because of repeatedly
downsampling and upsampling) in the merged maps, a 3$\times$3 covolution
filter is used to generate feature maps with the highest resolution and
accurate spatial information; therefore we can use them to do final
predictions.

##  4 Our Framework for Forest analysis

In this section, we explain how we applied Mask-RCNN to predict the location
of dead trees in forests. The framework structure can be referred to Figure 1.
All experiments are conducted using Tensorflow-GPU docker on a Mac Air with
eight NVIDIA Tesla V100 GPUs in Amazon Web Services (AWS).

###  4.1 Dataset collection and pre-processing

The dataset was collected by aerial photography on May 15th 2019, from the
Wood of Cree in Scotland, and was stored in a 10 GB merged tiff file. As the
file was too large for much of our software to process, we divided it into
40$\times$10 patches, each of 800$\times$800 pixels, using the Gridsplitter
plugin for the QGIS software. We then created our synthetic dataset from the
raw data: we extracted over 300 images of dead trees as foreground objects and
randomly resized, rotated and changed the lightness before they were placed
into 63 candidate backgrounds, also extracted from the raw data. Next we
produced a mask by filling random colours over the dead trees and black for
the background. We then retrieved $x$ and $y$ coordinates of the target
objects and stored them in a COCO format annotation file. Eventually, the
dataset was extended from 225 patches to 5000 patches and produced 5000
annotations in only a few hours.

Figure 9: Structure of Resnet-FPN, with the bottom-up pathway on the left, the
top-down pathway on the right. [37]. Image from [link] Figure 10: A raw image,
and the four corresponding feature maps. The lighter-coloured “hot zones”
indicate the locations of the features we’re interested in: shapes and colours
of dead trees.

###  4.2 How the backbone generates feature maps

Once the dataset is fed to our convolutional backbone, it produces four
feature maps for each picture in the dataset (see Figure 10).

The feature maps are then scanned by RPN in a sliding-window fashion, finding
anchors (pre-defined bounding boxes with fixed height and width) of high
objectness score, i.e. the probability that an anchor contains an object or a
non-object. Postive anchors have high objectness scores when they contain an
object.

Next, they are measured by value of 0.7 or more for the ”Intersection over
Union” (IoU) score. IoU is calculated as the ratio of the area of the overlap
(intersection) between the pre-defined bounding box and the ground truthbox
with the total area covered by both (union).

Note that because positive anchors do not always include the whole boundary of
their target object, they are examined by a regressor to refine the anchor
location. However, if several positive anchors overlap too much, we keep the
one with the highest objectness score and discard the rest. This method is
referred to as Non-max Suppression. Figure 11 shows the steps to merge the
anchor boxes.

Before the region proposals are sent to the second stage, the algorithm scans
and cuts out regions of interest (ROIs) from feature maps, as suggested by
region proposals. ROIs are then scaled to a pre-defined, reduced dimension
through the ROI pooling layer. (see Figure 6 )

Simultaneously, the selected ROIs are moved to the spatial locations
corresponding to the original picture in the ROI alignment layer, and then are
individually sent for pixel-wise classification in the FCN layer, resulting in
a prediction result as in Figure 12.

Figure 11: Step by step detection from finding postive anchor boxes to trim
down low confidence objects. Figure 12: Mask-RCNN prediction result.

###  4.3 Model fine-tuning

To make our training sessions faster and to more precisely predict target
objects, we made some modifications to the Mask-RCNN model.

First, we adopted a deeper model, i.e. training a ‘head + all + all’ model
instead of a single-head model, and we redefined the learning rate for each
training session. We set different learning rates based on optimizers, as in
[38]. Our choices of parameters are summarised in Table 1. Column with bold
type shows the a new, fine-tune value is assigned, compared to its last model.
Models are seperated to 4 groups with more similarity within group and are
allowed to be compared.

Model |  | Pre-
---
trained
net
| Acti-
---
vation
func-
tion
| Learn-
---
ing
rate
| Opti-
---
mizer
| Drop-
---
out
| Data-
---
set
| Fine-
---
tune-1
| Image-
---
Net
Relu |  |  $10^{-4}$,
---
$10^{-5}$,
$10^{-6}$
Adam | Yes | 2300
| Fine-
---
tune-2
| Image-
---
Net
| Leaky
---
Relu
|  $10^{-4}$,
---
$10^{-5}$,
$10^{-6}$
Adam | Yes | 2300
| Fine-
---
tune-3
| Image-
---
Net
| Leaky
---
Relu
|  $10^{-4}$,
---
$10^{-5}$,
$10^{-6}$
Adam | Yes | 3600
| Fine-
---
tune-4
| Image-
---
Net
Relu |  |  $10^{-4}$,
---
$10^{-5}$,
$10^{-6}$
Adam | Yes | 3600
| Fine-
---
tune-5
COCO | Relu |  |  $10^{-4}$,
---
$10^{-5}$,
$10^{-6}$
SDG | Yes | 3600
| Fine-
---
tune-6
COCO | Relu |  |  $10^{-4}$,
---
$10^{-5}$,
$10^{-6}$
SDG | No | 3600
| Fine-
---
tune-7
COCO | Relu |  |  $\mathbf{10^{-3}}$,
---
$\mathbf{10^{-3}}$,
$\mathbf{10^{-3}}$
SDG | No | 3600
| Fine-
---
tune-8
| Image-
---
Net
Relu |  |  $10^{-4}$,
---
$10^{-5}$,
$10^{-6}$
SDG | No | 4000

Table 1: Comparison of settings for our fine-tuning of Mask-RCNN.Column with
bold type shows the a new, fine-tune value is assigned, compared to its last
model. Models are seperated to four groups with more similarity within group.

Secondly, to avoid overfitting, we used data augmentation, i.e. the algorithm
accepts the original dataset, randomly transforms it and returns only the new,
transformed data [39][40]. We applied rotation by 0, 90, 180, and 270 degrees,
tuned brightness and contrast between 0.9 to 1.1, applied horizontal flipping
and clipped corners.

Thirdly, to make the loss function converge faster, we changed our activation
function from Relu to Leaky Relu. The latter avoids the problem with “dying
neurons” which Relu activation functions suffer from: when a neuron “dies”,
whichever input it receives, its output will be equal to 0; in that case a
neuron loses prediction power and stops learning. Leaky Relu, however, has a
small but non-zero gradient for inputs below 0.

Fourthly, we applied a dropout function in the FPN layer to randomly drop out
input nodes before doing classification and predicting bounding boxes. This
decreases random correlations within neurons. It serves a similar purpose as
data augmentation, i.e. to prevent overfitting [41].

Fifthly, we changed the optimizer from stochastic gradient descent (SGD) to
Adaptive Moment Estimation (Adam)[42]. An optimizer aims to update parameters
$\theta$ in the negative gradient direction to minimize the loss. Studies show
that Adam has less generalization error compared to SGD as well as being able
to tolerate a noisy or sparse gradient. The algorithm is denoted

$\displaystyle m_{t}$ $\displaystyle=\beta_{1}m_{t-1}+(1-\beta_{1})g_{t}$
$\displaystyle v_{t}$ $\displaystyle=\beta_{2}v_{t-1}+(1-\beta_{2})g_{t}^{2}$

where $m_{t}$ and $v_{t}$ are the decaying averages of the past gradient and
past square gradient respectively, and the default value of each of the two
$\beta$s is close to 1. We then update the estimates of both to deal with the
bias caused by the initial averages $v=0$ and $m=0$:

$\displaystyle\hat{m}_{t}$ $\displaystyle=\frac{m_{t}}{1-\beta^{t}_{1}}$
$\displaystyle\hat{v}_{t}$ $\displaystyle=\frac{v_{t}}{1-\beta^{t}_{2}}.$

Because we initialize averages with zeros, the estimators are biased towards
zero. However, the author of Adam optimizer used the above formulas to correct
the estimators. [42]

Then, using the updated parameters, the Adam update rule is:

$\theta_{t+1}=\theta_{t}-\frac{\eta}{\sqrt{\hat{v}_{t}}+\epsilon}\hat{m}_{t},$

where $\eta$ is the learning rate, and we choose (as part of our tuning of the
model) $\epsilon=10^{-8}$.

Finally, we compare prediction results from fine-tune experiments (ex,
activation function, dropout function, optimizer) to model without fine tune
by visualizing loss history against epoch and a table reveals mean average
precison (mAP) [43].

###  4.4 Loss function of Mask-RCNN

Loss functions are used to evaluate how well a model performs [18] \- in our
case we aim to minimise a loss function by optimising weights for three
different tasks: classification, bounding box prediction and mask prediction.
We use a loss function which is the sum of a loss function for each of these
tasks:

$L=L_{\mathrm{cls}}+L_{\mathrm{box}}+L_{\mathrm{mask}}$

where $L_{\text{cls}}$ is the same as in Faster RCNN:

$L_{\mathrm{cls}}(P_{i},P^{*}_{i})=-\frac{1}{N_{\mathrm{cls}}}\left[P^{*}_{i}\log(P_{i})+(1-P^{*}_{i})\log(1-P_{i})\right]$

where $P_{i}$ is the _predicted_ probability that the $i$th anchor is a target
object. In Matterport Mask RCNN, the total number of anchors is 256 (positive
and negative). The $i$th anchor above “$i$th” means a random number from 1 to
256.

If an anchor contains target objects, the _ground truth label_ $P^{*}_{i}$ is
equal to 1, and otherwise it is equal to 0. In addition, the function is
normalised by $N_{\mathrm{cls}}$, i.e. the mini-batch size. A batch is the
number of training samples in a given iteration. Overall, the classification
loss $L_{\mathrm{cls}}$ is a log-loss over two classes (i.e. target object and
not target object) which can be thought of as a binary classification loss
(cross-entropy loss). In addition, $L_{\mathrm{box}}$ is also as in Faster
RCNN:

$L_{\mathrm{box}}=\frac{\lambda}{N_{\mathrm{cls}}}P_{i}\,R((t_{i}-t^{*}_{i})_{i\in\\{\mathrm{x,y,w,h}\\})}$

where $R$ is a robust $L_{1}$ loss:

$R(\mathbf{t}^{u}-\mathbf{v})=\sum_{i\in\\{x,y,w,h\\}}\mathrm{smooth}_{L_{1}}(t^{u}_{i}-v_{i})$

and $t_{i}^{u}$ for $i\in\\{x,y,w,h\\}$ are the position and dimension
parameters of the predicted bounding box, $v_{i}$ the parameters of the ground
truth box, $\lambda$ is a balancing parameter which is set to be around 10 in
the paper (so that both losses $L_{\mathrm{cls}}$ and $L_{\mathrm{box}}$ are
roughly equally weighted), and

$\mathrm{smooth}_{L_{1}}(s)=\left\\{\begin{matrix}0.5s^{2}&\text{ if
}|s|<1,\\\ |s|-0.5&\text{ otherwise.}\end{matrix}\right.$

We also have parameterised coordinates $t_{i}$ and $t^{*}_{i}$,

$\displaystyle(t_{i})$ $\displaystyle=(t_{x},t_{y},t_{w},t_{h})$
$\displaystyle(t^{*}_{i})$
$\displaystyle=(t^{*}_{x},t^{*}_{y},t^{*}_{w},t^{*}_{h})$

where each individual parameterised coordinate is calculated given the
predicted bounding box, ground truth box and anchor box coordinates as

$\begin{matrix}t_{x}=(x-x_{\mathrm{a}})/w_{\mathrm{a}},&t_{y}=(y-y_{\mathrm{a}})/h_{\mathrm{a}},\\\
\vspace{0.2cm}t_{w}=\log(w/w_{\mathrm{a}}),&t_{h}=\log(h/h_{\mathrm{a}})\\\
t_{x}^{*}=(x^{*}-x_{\mathrm{a}})/w_{\mathrm{a}},&t_{y}^{*}=(y^{*}-y_{\mathrm{a}})/h_{\mathrm{a}},\\\
t_{w}^{*}=\log(w^{*}/w_{\mathrm{a}}),&t_{h}^{*}=\log(h^{*}/h_{\mathrm{a}}),\end{matrix}$

where $(x_{\mathrm{a}},y_{\mathrm{a}},w_{\mathrm{a}},h_{\mathrm{a}})$ are the
coordinates of the anchor box.

Finally, $L_{\mathrm{mask}}$ is the average binary cross-entropy loss which is
calculated when the $k$th mask is the associated ground truth mask:

$L_{\mathrm{mask}}=-\frac{1}{m^{2}}\sum_{1\leq i,j\leq m}\left[y_{ij}\log
y_{ij}^{k}+(1-y_{ij})\log(1-y_{ij}^{k})\right]$

where $m\times m$ is the mask dimension for each ROI in each class, $y_{ij}$
is the ground truth label (0 or 1) in cell $(i,j)$ , while $y^{k}_{ij}$ is the
predicted probability that the same cell belongs to the $k$th class [18].

##  5 Experiments

###  5.1 Implement on Cloud Server with GPUs

A cloud server provides a virtual environment to perform computation which can
be accessed remotely via internet. The major benefit is to overcome hardware
limitation, for example, lack of disk space to store data or computing power
to train models. In this study, in order to deal with the enormous amount of
data and complicated algorithms, a cloud server(AWS) is chosen to provide the
environment needed for the AI program (Mask RCNN).

Once an AWS computing instance has been registered, a private key file is
given access to cloud server. Then, a Tensorflow-GPU docker file is used to
activate eight GPUs as well as giving the user a temporary environment with
numerous pre-installed python libraries. After the training session is done,
the post-trained files are transfered back to cloud server before closing
docker environment and retrieve files to local computer.

The training process takes about 4-6 hours for a dataset of 2000 to 4000
images.

###  5.2 Evaluation of prediction performance

As well as the loss function, we use COCO mean average precision (COCO mAP) to
evaluate prediction performance [43]. This indicator is widely used in
evaluating models for object detection. Firstly, after IoU decides the ROIs,
we calculate @mAP based on the prediction performance of those ROIs (ignoring
objects with an IoU score lower than 50%). Then we calculate the mean
precision, recall and F1. Then the AP is calculated as the area under the
precision-recall curve: the so-called Area Under Curve (AUC). The curve shows
the trade-off relationship between precision and recall. When the area is
larger, the predicted results are more accurate (numbers of false positives
and false negatives are lower) and more positive cases are detected (number of
false nagative is lower).

In addition, in COCO AP, the AP averaged over IoU thresholds is calculated,
i.e. AP50, AP75 and AP95. In COCO mAP, an average for a 101-point interpolated
AP is calculated. This means recalls are divided into 101 points and their
maximized precisions are averaged:

$\displaystyle\mathrm{AP}$
$\displaystyle=\frac{1}{101}\sum_{r\in\\{0,0,\cdots,1,0\\}}\mathrm{AP}_{r}$
$\displaystyle=\frac{1}{101}\sum_{r\in\\{0,0,\cdots,1,0\\}}P_{\mathrm{interp}}(r).$

To make an even more direct indicator, we also calculate the ROI number of two
types of target objects (Figure 13) which will be revealed in the experimental
results.

###  5.3 Experimental results and analysis

Because CNN can easily suffer from overfitting, we started training our model
with a dataset of 2 000 images and a validation set of 300 images, to which we
applied many transformations. We then looked to increase the size of the
dataset gradually to find the training set which gives us the most general
model. We trained each model with 100 epochs with different hyperparameter
settings.

For the first pair of model, Fine-tune-1 and Fine-tune-2, the results showed
that for the two types of targets (types A and B, see Figure 13) present in
the data, ROI detection did not perform well for predicting type B objects
(round shape) in the model “Fine-tune-1”. However, the results also showed the
opposite trend in model “Fine-tune-2”, i.e. type B objects were detected
reliably but few type A objects were. To resolve this, we decided to include
more data in the next trial.

Figure 13: The model “Fine-tune-1” can reveal the type A pattern (left) but
not type B (right), and “Fine-tune-2” can reveal the type B pattern but not
type A.

For the next models, we reveal the evaluation statistics in Table 2. Note that
we make an abbreviation of ”fine-tune” to ”FT”, for example, “FT3” is the
model “Fine-tune-3”. In addition, we indicate the best three performances for
each measurement in bold:

Model | FT3 | FT4 | FT5 | FT6 | FT7 | FT8
---|---|---|---|---|---|---
AP50 | 0.39 | 0.57 | 0.43 | 0.54 | 0.60 | 0.38
AP75 | 0.14 | 0.53 | 0.13 | 0.48 | 0.56 | 0.10
mPrecision | 0.62 | 0.78 | 0.62 | 0.76 | 0.82 | 0.74
mRecall | 0.35 | 0.39 | 0.35 | 0.40 | 0.41 | 0.34
mF1 | 0.29 | 0.32 | 0.28 | 0.33 | 0.34 | 0.28
Type A (star) | 40 | 59 | 33 | 69 | 32 | 53
Type B (round) | 71 | 24 | 64 | 74 | 24 | 64

Table 2: Comparison of the performance of the models. The best three
performances for each measurement are in bold. The final two rows show
bounding box numbers for the two types of target objects.

The table shows that results of FT6 are the best among all of the fine-tuned
models. In addition, in Figure 14 we plotted the loss results to see if FT6
had the lowest training and validation loss.

Figure 14: Training and validation loss for all models.

Because FT7 has the lowest training loss and FT6 has lowest validation loss in
Figure 14, we compared FT6 with FT7.

Figure 15: Total loss for the models FT6 and FT7.

FT7 had a constant higher learning rate of $10^{-3}$ which is the closest of
all our models to Matterport’s Mask RCNN default settings [31], while FT6
followed the fast-medium-slow learning rate of $10^{-4}$-$10^{-5}$-$10^{-6}$.
The results in Figure 15 show that FT7 suffered from both fluctuating loss and
overfitting. Overfitting happens when the training loss is less than the
validation loss which means the algorithm learns local features appearing in
the training set but not the validation set. This indicates that a fast and
fixed learning rate in FT7 has no advantage over a learning rate which follows
a convex pattern.

In addition, there were two major trends in Figure 14: a big shift in the
training loss, and two peaks in the validation loss. Because all models had
different settings, especially in the different groups, we assumed that the
enlarged data set helped decrease the training loss. In addition, for the
validation loss, two peaks were made by FT1 and FT4, and so we believe that
this may be because of their model setting (Relu + Adam + Dropout).

The distributions of precision, recall and F1 from the model FT6 are shown
below in Figure 16. The diagrams show that FT6 has very high precision but a
rather low recall. This means number of false negatives is higher than the
number of false postive.

Figure 16:  Distributions of precision, recall and F1 score for FT6
(respectively). The x-axis counts the number of predicted bounding boxes while
the y-axis records the precision, recall and F1 scores respectively.

##  6 Discussion

###  6.1 Advantages and disadvantages

Advantages of our method

  1. 1.

Efficently produces images and annotations:

The synthetic method was designed by Adam Kelly to avoid a time-consuming
image collection and annotation process. In Adam’s experiment [9], he created
2000 images containing cigarette butts and their annotations in 20 minutes
from only 25 samples. Because we had small numbers of data and no longer
collected dead tree images when we started to train our algorithm, we chose
the synthetic method rather than trying to find more remote sensor data.

  2. 2.

Less complex to apply compared to MaxEnt algorithm in CanHeMon project:

According to the image pre-processing section of the CanHeMon project, there
are 14 potential predictor variables related to texture layer, containing
significant image information. Then, after a series of background and AUC
setting, the MaxEnt algorithm and principal component analysis are applied to
reduce the number of superfluous variables. However, when applying CNN based
algorithm, the variables used as inputs and ouputs are rather simple: raw
images and their labels. There is no manual elimination process for the
variables/features but the algorithm automatically adjusts feature weights
after backpropagation in each iteration.

  3. 3.

False positive cases are rare:

The CanHeMon project suffered from a high rate of false positives, which does
not occur in our Mask RCNN-based model. In the CanHeMon experiment,
researchers manually removed false postive cases and fed the data back to the
training session. However, the mask results from our model rarely showed false
postive cases, although there were still many areas of dead trees which were
not detected.

We found that a decrease in the IoU threshold can help detect more dead trees,
and therefore decrease the number of false negatives.

Changing a hyperparameter in this way requires less manual effort than
manually dropping FP cases.

Disadvantages of our method

  1. 1.

Data scarcity:

We gathered all of the images of dead trees we could from both online
resources and our partner company 2Excel. However, we did not generate enough
data to overcome the problems caused by data scarcity, as dead trees are an
uncommon object, whereas most image datasets collect everyday objects (cars,
people, etc.). The main problem was that there were several shapes of dead
tree in our raw data (see the two types in Figure 13) but we could not
generate a sufficient number of images to deal with this variety.

  2. 2.

The process of producing COCO format annotations did not give clear boundaries
in synthetic images:

The synthetic method we used sometimes led to target objects being divided
into more than one piece by another object overlapping them. For example, see
the red region in Figure 17. We had to manually drop the label when this
happened, which is time consuming, reducing the benefit in time savings we
gained from using the synthetic method. The overlapping type B dead trees give
us inaccurate bounding box predictions.

Figure 17: A target object in this synthetic image is divided in two by
another object placed on top of it.

  3. 3.

Predicted bounding boxes failing to accurately cover a whole region of dead
wood

Our synthetic method can make target objects overlap, and if they belong to
type A (star shaped) then the algorithm will fail to distinguish the objects
from each other, as in Figure 18. At best, they are bound to a larger group
and given a single label.

Figure 18: Several target objects of type A being assigned to a single
bounding box.

  4. 4.

Reliability of mAP and other statistics

The mAP scores in our results were still not as reliable as the ones produced
and evaluated using a raw (not synthetic) validation set. The synthetic
validation set was necessary because of the scarce data, but it would be ideal
to use only raw images for validation if possible.

###  6.2 Suggestions for future research

In this research, we only trained our model with the company’s aerial dataset
instead of the Landsat dataset; however, we welcome to reseachers apply Mask
RCNN to the Landset datset to make comparison with our results. In addition,
we did not take into account other variables such as the forest phenology in
our model, because the aerial dataset was collected in summer time (May 15th
2019) when living trees are easily distinguished from dead trees. However, we
cannot exclude the possibility that Mask RCNN detects and labels as “dead”
rare species of trees who naturally have white leaves. Because of that, we
encourage researchers to apply other possible features in their models rather
than only relying on an image processing algorithm.

Secondly, to mitigate the disadvantage of data scarcity, we enlarged our
dataset and produced labels automatically using a synthetic method. However,
we pointed out some issues with this method in subsection 6.1, although recent
updates to the code we were using have dealt with some of the issues we saw in
item 2 of the disadvantages in subsection 6.1. As Adam Kelly suggested [9], it
is still necessary to evaluate the advantages and disadvantages of using
synthetic data to train a model.

The issue with bounding box inaccuracy is the major contributor to the high
loss in RPN. A possible solution is using a single FCN to classify each pixel
if the instance-level segmentation is not necessary. Then, we could change the
loss function from bbox + mask to FCN loss, decreasing the probability of high
loss and enabling a more flexible and accurate prediction because mask
prediction is not limited by the selected ROIs’ areas.

One other notable issue we experienced was how to increase the reliability of
statistics (accuracy, mAP, mRecall and mF1) which were produced using the
synthetic validations set to be consistent with the model’s predictive power
on raw images. However, if we were able to collect a sufficiently large,
manually annotated dataset, the statistics would be much more reliable. In
addition, if the synthetic training set is replaced by manually-annotated raw
data, it could increase the predictive power since the synthetic data may not
include features which are present in the raw data.

In the conventional training process, a single optimizer is used in a model.
However, to take advantage of the merits of different optimizers, the author
in [44] used ADAM first, and when the loss stopped decreasing, replaced ADAM
with SGD. The results obtained by doing this were better than those obtained
using either optimizer on its own. After reading this, we performed the same
experiment with our model. As Figure 19 shows, replacing SGD with ADAM at the
30th epoch led to a sudden drop in the loss function. As it is hard to compare
directly with the other models we tested, we did not include this model in
Table 1, but the results of our experiment are encouraging, so we would
suggest that future research should look into using multiple optimizers for
training.

Figure 19: The loss functions obtained by using an SGD optimiser until epoch
30, and then using ADAM.

Figure 20 shows more of our test results. Mostly the detection framework works
well on these images while the training process was carried out on the same
dataset. However, there is yet a gap between the expected detection accuracy
and the performance of our deep learning framework, and further improvement is
needed to extend our framework to practical applications.

Figure 20: More test results on aerial images.

The limits can be ascribed to two aspects. First, deep neural networks has
been reported not robust to noises and other similar perturbations (adversary
attacks)[45], while DNNs are somehow lack of explicit explanability unlike
typical statistic approaches that model the data space explicitly
[46][47][48][49]. Secondly, DNNs are date-driven and heavily rely on training
datasets [50][51][52][53]. In our experiments, we have a limited amount of
data samples. A typical remedy to this issue can be based on the neural
network models of Generative Adversarial Networks (GANs)[54], which can
generate more similar images similar to real ones to improve the training
process.

Targeting at dead wood detection, our work is based on the assumption that the
dead tree dataset collected in the spring-summer period have different colors
in leaves. However, this may vary according to tree species and some rare
specious of trees may naturally have white leaves during four seasons. A
future work can include forest phenology and classify the species of trees
before the dead wood detection.

##  7 Conclusion

In conclusion, we present a new framework for automated dead tree detection
from aerial images using a re-trained Mask RCNN approach, with a transfer
learning scheme. We apply our framework to our aerial imagery datasets, and
compare eight fine-tuned models. The mean average precision score (mAP) for
the best of these models reaches 54%. We also are able to automatically
produce mask visualizations to label the dead trees in an image, so that the
number of dead trees in a region can be automatically counted. Such aerial
image-based forest analysis can give an early diagnosis of the conditions of
forests and identify potential risk factors, and provide an potential solution
to avoid of climate change caused disasters such as the recent forest fires in
the California and Australia.

## Acknowledgements

This work was supported in part by the EPSRC grant (EP/P009727/1).

## References

  * [1] “Ancient tree guide 5: Trees and climate change,” _WoodlandTrust_ , 2008.
  * [2] S. Florence, “Half a billion euros to flow into germany’s forests,” _EURACTIV_ , 2019, Aug 30.
  * [3] A. Mike, “Self-induced disaster: California fires the direct result ofshortsighted environmentalist policies that prohibit forest management,” _ECOLOGY NEWS_ , 2019, Oct 11.
  * [4] S. Joe and M. Bruno, “This summer uncontrolled fires have raged acrossthe uk and europe as never before. what, or who, sparked them? andmight we discover the danger is far from over when the flames aredoused?” _The Telegraph_ , 2019, Aug 3.
  * [5] P. S. Beck, L. Martínez, M. d. L. Sanchez, Y. Chemin, G. Caudullo _et al._ , “The canopy health monitoring (canhemon) project,” 2019.
  * [6] L. T and J. S, “Transfer learning,” in _Handbook of research on machine learning applications and trends: algorithms, methods, and techniques_. IGI Global, 2010, pp. 242–264.
  * [7] T.-Y. Lin, M. Maire, S. Belongie, J. Hays, P. Perona, D. Ramanan, P. Dollár, and C. L. Zitnick, “Microsoft coco: Common objects in context,” in _European conference on computer vision_. Springer, 2014, pp. 740–755.
  * [8] Y. Guo, Y. Liu, T. Georgiou, and M. S. Lew, “A review of semantic segmentation using deep neural networks,” _International journal of multimedia information retrieval_ , vol. 7, no. 2, pp. 87–93, 2018.
  * [9] “Complete guide to creating coco datasets.” [Online]. Available: https://github.com/akTwelve/cocosynth
  * [10] D. S. Boyd and F. M. Danson, “Satellite remote sensing of forest resources: three decades of research development,” _Progress in Physical Geography: Earth and Environment_ , vol. 29, no. 1, pp. 1–26, 2005.
  * [11] K. S. Yankovich, E. P. Yankovich, and N. V. Baranovskiy, “Classification of vegetation to estimate forest fire danger using landsat 8 images: Case study,” vol. 2019, pp. 1–14, 2019.
  * [12] J. Miettinen, H.-J. Stibig, and F. Achard, “Remote sensing of forest degradation in southeast asia—aiming for a regional view through 5–30 m satellite data,” vol. 2, pp. 24–36, 2014.
  * [13] P. Borrelli, L. A. S. Rondón, and B. Schütt, “The use of landsat imagery to assess large-scale forest cover changes in space and time, minimizing false-positive changes,” vol. 41, pp. 147–157, 2013.
  * [14] “Landsat data access.” [Online]. Available: https://www.usgs.gov/land-resources/nli/landsat
  * [15] “Landsat—earth observation satellites,” _U.S. Geological Survey (USGS)_ , 2016. [Online]. Available: https://pubs.er.usgs.gov/publication/fs20153081
  * [16] R. Girshick, J. Donahue, T. Darrell, and J. Malik, “Rich feature hierarchies for accurate object detection and semantic segmentation.”
  * [17] R. Girshick, “Fast r-cnn,” in _Proceedings of the 2015 IEEE International Conference on Computer Vision (ICCV)_ , ser. ICCV ’15. USA: IEEE Computer Society, 2015, p. 1440–1448. [Online]. Available: https://doi.org/10.1109/ICCV.2015.169
  * [18] S. Ren, K. He, R. Girshick, and J. Sun, “Faster r-cnn: Towards real-time object detection with region proposal networks,” in _Advances in neural information processing systems_ , 2015, pp. 91–99.
  * [19] J. Redmon, S. Divvala, R. Girshick, and A. Farhadi, “You only look once: Unified, real-time object detection.”
  * [20] W. Liu, D. Anguelov, D. Erhan, C. Szegedy, S. Reed, C.-Y. Fu, and A. C. Berg, “Ssd: Single shot multibox detector,” pp. 21–37, 2016.
  * [21] B. Zoph, V. Vasudevan, J. Shlens, and Q. V. Le, “Learning transferable architectures for scalable image recognition,” 2018.
  * [22] Z.-Q. Zhao, P. Zheng, S. tao Xu, and X. Wu, “Object detection with deep learning: A review.”
  * [23] J. Huang, V. Rathod, C. Sun, M. Zhu, A. Korattikara, A. Fathi, I. Fischer, Z. Wojna, Y. Song, S. Guadarrama, and K. Murphy, “Speed/accuracy trade-offs for modern convolutional object detectors.”
  * [24] K. He, G. Gkioxari, P. Dollár, and R. B. Girshick, “Mask R-CNN,” _CoRR_ , vol. abs/1703.06870, 2017.
  * [25] J.-D. Sylvain, G. Drolet, and N. Brown, “Mapping dead forest cover using a deep convolutional neural network and digital aerial photography,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 156, pp. 14 – 26, 2019.
  * [26] S. Jiang, W. Yao, and M. Heurich, “Dead wood detection based on semantic segmentation of vhr aerial cir imagery using optimized fcn-densenet,” vol. XLII-2/W16, pp. 127–133.
  * [27] T. Kattenborn, J. Eichel, and F. E. Fassnacht, “Convolutional neural networks enable efficient, accurate and fine-grained segmentation of plant species and communities from high-resolution uav imagery,” vol. 9, 2019.
  * [28] A. Mazza, F. Sica, P. Rizzoli, and G. Scarpa, “Tandem-x forest mapping using convolutional neural networks,” vol. 11, p. 2980.
  * [29] Y. Wei and L. Tao, “Efficient histogram-based sliding window,” 2010.
  * [30] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 1097–1105.
  * [31] “Matterport mask r-cnn.” [Online]. Available: https://github.com/matterport/Mask_RCNN
  * [32] A. Albert, J. Kaur, and M. C. Gonzalez, “Using convolutional networks and satellite imagery to identify patterns in urban environments at a large scale,” in _Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining_ , 2017, pp. 1357–1366.
  * [33] B. Tan, Y. Zhang, S. J. Pan, and Q. Yang, “Distant domain transfer learning,” in _Thirty-First AAAI Conference on Artificial Intelligence_ , 2017.
  * [34] R. Yamashita, M. Nishio, R. K. G. Do, and K. Togashi, “Convolutional neural networks: an overview and application in radiology,” _Insights into imaging_ , vol. 9, no. 4, pp. 611–629, 2018.
  * [35] S. Dey, _Hands-On Image Processing with Python: Expert techniques for advanced image analysis and effective interpretation of image data_. Packt Publishing, 2018.
  * [36] V. Dumoulin and F. Visin, “A guide to convolution arithmetic for deep learning,” _arXiv preprint arXiv:1603.07285_ , 2016.
  * [37] T.-Y. Lin, P. Dollár, R. Girshick, K. He, B. Hariharan, and S. Belongie, “Feature pyramid networks for object detection,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 2117–2125.
  * [38] Y. Wu, L. Liu, J. Bae, K.-H. Chow, A. Iyengar, C. Pu, W. Wei, L. Yu, and Q. Zhang, “Demystifying learning rate polices for high accuracy training of deep neural networks,” _arXiv preprint arXiv:1908.06477_ , 2019.
  * [39] “Image preprocessing: Keras imagedatagenerator class.” [Online]. Available: https://keras.io/preprocessing/image/
  * [40] K. Li, G. Cheng, S. Bu, and X. You, “Rotation-insensitive and context-augmented object detection in remote sensing images,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 56, no. 4, pp. 2337–2348, 2017.
  * [41] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, “Dropout: a simple way to prevent neural networks from overfitting,” _The journal of machine learning research_ , vol. 15, no. 1, pp. 1929–1958, 2014.
  * [42] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [43] K. Oksuz, B. Can Cam, E. Akbas, and S. Kalkan, “Localization recall precision (lrp): A new performance metric for object detection,” in _Proceedings of the European Conference on Computer Vision (ECCV)_ , 2018, pp. 504–519.
  * [44] N. S. Keskar and R. Socher, “Improving generalization performance by switching from adam to sgd,” _arXiv preprint arXiv:1712.07628_ , 2017.
  * [45] R. Jiang and D. Crookes, “Shallow unorganized neural networks using smart neuron model for visual perception,” _IEEE Access_ , vol. 7, pp. 152 701–152 714, 2019.
  * [46] R. Jiang, A. T. Ho, I. Cheheb, N. Al-Maadeed, S. Al-Maadeed, and A. Bouridane, “Emotion recognition from scrambled facial images via many graph embedding,” _Pattern Recognition_ , vol. 67, pp. 245–251, 2017.
  * [47] R. M. Jiang, D. Crookes, and N. Luo, “Face recognition in global harmonic subspace,” _IEEE Transactions on Information Forensics and Security_ , vol. 5, no. 3, pp. 416–424, 2010.
  * [48] R. M. Jiang, D. Crookes, N. Luo, and M. W. Davidson, “Live-cell tracking using sift features in dic microscopic videos,” _IEEE Transactions on Biomedical Engineering_ , vol. 57, no. 9, pp. 2219–2228, 2010.
  * [49] R. M. Jiang, A. H. Sadka, and D. Crookes, “Multimodal biometric human recognition for perceptual human–computer interaction,” _IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews)_ , vol. 40, no. 6, pp. 676–681, 2010.
  * [50] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” _nature_ , vol. 521, no. 7553, pp. 436–444, 2015.
  * [51] G. Storey, R. Jiang, and A. Bouridane, “Role for 2d image generated 3d face models in the rehabilitation of facial palsy,” _Healthcare technology letters_ , vol. 4, no. 4, pp. 145–148, 2017.
  * [52] Z. Jiang, P. L. Chazot, M. E. Celebi, D. Crookes, and R. Jiang, “Social behavioral phenotyping of drosophila with a 2d–3d hybrid cnn framework,” _IEEE Access_ , vol. 7, pp. 67 972–67 982, 2019.
  * [53] G. Storey, R. Jiang, S. Keogh, A. Bouridane, and C.-T. Li, “3dpalsynet: A facial palsy grading and motion recognition framework using fully 3d convolutional neural networks,” _IEEE Access_ , vol. 7, pp. 121 655–121 664, 2019.
  * [54] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” in _Advances in neural information processing systems_ , 2014, pp. 2672–2680.


# Non-exponential decay in classical stochastic processes with memory††thanks:
e-mails<EMAIL_ADDRESS><EMAIL_ADDRESS>
Maciej Rybczyński and Zbigniew Włodarczyk  Institute of Physics, Jan
Kochanowski University, PL-25406 Kielce, Poland

###### Abstract

The initial time-dependence of a state in circumstances where it makes
transitions to, or decay to, a second state has been investigated. In
classical stochastic processes, the observed time dependence of transition or
decay proportional to $t^{2}$ is attributed to the noise with memory. In
contrast to quantum mechanics, the quadratic form of initial decay is unable
to decelerate the evolution of the system.

##  1 Introduction

In many stochastic processes the survival probability $S(t)$ attracts
considerable attention [1, 2, 3, 4, 5]. The standard exponential decay law
$S(t)=e^{-t/\tau}$ is never exact. For short times, $t\ll\tau$, the decay law
is typically quadratic

$S\left(t\right)\simeq 1-\left(\frac{t}{\tau}\right)^{2}.$ (1)

Obviously a central part of such analysis is the use of the $t^{2}$ factor
(or, at last in principle, $t^{\alpha}$ with $\alpha>1$) to produce a slowing
down of decay by measurement. As a consequence, for frequently repeated
measurements at very small time intervals, the Zeno effect, that is the
freezing of state in its unstable initial configuration, takes place [6, 7,
8]. If $S\left(t\right)$ oscillates around exponential decay law, also the
inverse Zeno effect, that is an increased decay rate by measurements, is
possible [9]. Finally, at large times $t\gg\alpha\tau$ a power law sets in,
$S(t)\simeq t^{-\alpha}$, where $\alpha>0$ depends on the details of the
interaction that rules that particular decay [10] 111Non-exponential decay of
unstable systems can also be understood and described in a most economical and
adequate way in terms of non-extensive statistics [11]. It is the fact that
$t^{2}-$ dependence for the quantum case leads to the quantum Zeno prediction,
and should be regarded as the quantum mechanical norm. There is, thus, a
fundamental mismatch between quantum physics and classical physics, both in
the very specific point of the $t^{2}-$ dependence versus the $t-$ dependence
for short time, and in the more general point that decay processes are
exponential in classical case, but not in the quantum one. In this work we
discuss classical processes influenced by noise with memory, which show
$t^{2}$-form of initial decay, and we show that the Zeno effect vanishes in
this case.

In Section 2 we recall the standard analysis of quantum decay. In section 3 we
construct master equation with stochastic noise which solution for small
enough $t$ show $t^{2}-$ form. Then, we discuss obtained solution (section 4),
and finally we advocate in section 5 that memory destroy independence and in a
classical system an unstable particle observed continuously decay quasi-
exponentially (the Zeno paradox does not exist). Section 6 summarizes and
concludes our work.

##  2 Quantum Zeno paradox

Simple derivation of the quantum Zeno effect is possible by considering time
behaviour of the state vector [12]. Let the state vector at time $t$ is
$e^{-iHt}|\Phi\rangle$, where $|\Phi\rangle$ is the state vector at time $t=0$
and $H$ is the Hamiltonian in units where the reduced Planck’s constant is set
to unity, $\hbar=1$. For $t$ small enough, it is possible to make a power
series expansion: $e^{-iHt}\simeq 1-iHt-H^{2}t^{2}/2+...$. The survival
probability is

$S(t)=|\langle\Phi|e^{-iHt}|\Phi\rangle|^{2}\simeq 1-(\Delta H)^{2}t^{2},$ (2)

where

$(\Delta H)^{2}=\langle\Phi|H^{2}|\Phi\rangle-\langle\Phi|H|\Phi\rangle^{2}$
(3)

is the variance of Hamiltonian.

Interrupted the interval $[0,t]$ by $n$ independent sequential measurements at
times $t/n$, the survival probability is

$S(t)\simeq\left[1-(\Delta H)^{2}\left(\frac{t}{n}\right)^{2}\right]^{n}$ (4)

and approaches 1 for $n\rightarrow\infty$. It illustrates the situation that
an unstable particle, if observed continuously, will never decay. Zeno effect
is a result of (a) a $t^{2}$-form of initial decay, which is the normal
quantum case (cf. Eq. (2)), and (b) the independence of sequential
measurements (cf. Eq. (4)) due to the fact, that in quantum mechanics,
measurements generally change the state of the system being measured [13].

##  3 Non-local master equation with memory

In classical case we consider time-dependence of a state in circumstance where
it may make transition to, or decay to, a second state. The rate of change of
the non-decay probability $P(t)$ at time $t$ and dependence on its history
(starting at $t=0$) is given by the non-local master equation

$\frac{dP(t)}{dt}=-\int_{0}^{t}\kappa(t-s)P(s)\,ds+\eta(t),$ (5)

where $\eta$ is the stochastic noise 222This approach to the evolution of
statistical system is based on the Mori-Zwanzig formalism in which memory
effects are taken into account through the introduction of the memory kernel
$\kappa\left(t\right)$ [14]. This means that the rate of changes of the state
$P\left(t\right)$ at time $t$ depends on its history (starting at $t=0$).
Recently, the Langevin equation in approach given by Eq. (5) was considered in
description of particles motion [15]. The memory effects are taken into
account through the introduction of the memory kernel

$\kappa(t,s)=\langle\eta(t)\eta(s)\rangle=\gamma f(t-s)$ (6)

with $\gamma=(\Delta\eta)^{2}=\langle\eta^{2}\rangle-\langle\eta\rangle^{2}$
being variance of $\eta$ (notice that $\langle\eta\rangle=0$). We assume that
$\eta$ is a Lorentzian noise with spectrum

$G(\omega)=\langle\eta^{2}\rangle\frac{\tau}{1+\omega^{2}\tau^{2}}$ (7)

with correlators given by
333$\langle\eta(0)\eta(t)\rangle=\int_{0}^{\infty}{G(\omega)e^{i\omega
t}d\omega}=\langle\eta^{2}\rangle e^{|t|/\tau}$

$f(t)=e^{-|t|/\tau}$ (8)

with $\tau$ being the memory time since it sets the time scale over which time
correlations of noise decay 444For $f(t-s)=\delta(t-s)$ we have standard
exponential decay law, $S(t)=e^{-t/\tau}$ [16].

The non-decay probability $P\left(t\right)$ is by itself a random variable
depending on the field fluctuations [17]. We define survival probability as
the expectation value of the non-decay probability. Averaging probability over
ensemble, for survival probability $S(t)=\langle P(t)\rangle$ we have

$\frac{dS(t)}{dt}=-\int_{0}^{t}\gamma f(t-s)S(s)\,ds.$ (9)

Differentiating both side of Eq. (9) we have

$\frac{d^{2}S(t)}{dt^{2}}+\frac{1}{\tau}\frac{dS(t)}{dt}+\gamma S(t)=0$ (10)

which, with two initial conditions:

$S\left(t=0\right)=0\mbox{~{}and~{}}\frac{dS\left(t\right)}{dt}\Biggl{|}_{t=0}=0,$

allow simply to solve the problem analytically. Survival probability is given
by

$S(t)=\frac{e^{-t(1+A)/(2\tau)}(-1+A)}{2A}+\frac{e^{-t(1-A)/(2\tau)}(1+A)}{2A}$
(11)

with $A=\sqrt{1-4\gamma\tau^{2}}$ (for formal solution of Eq. 5 see Appendix
A). Time dependence of the survival probability $S(t)$ for different
parameters $A$ is shown in Fig. (1). Only for large $t$ we observe exponential
behaviour, $S(t)\sim\exp[-t(1-A)/(2\tau)]$. For $t$ small enough we have

$S(t)\simeq 1-\frac{1}{2}\gamma t^{2}=1-\frac{1}{2}(\Delta\eta)^{2}t^{2}$ (12)

with the same behaviour as $S(t)$ given by Eq. (2). Comparison of parabolic
dependence given by Eq. (12) with initial time dependence of $S\left(t\right)$
is shown in Fig. (2). The variance of noise, characterized by
$(\Delta\eta)^{2}/2$ term, corresponds to $(\Delta H)^{2}$, and survival
probability $S(t)$ the same exhibit $1-t^{2}$ form [17].

Figure 1: Survival probability $S(t)$ given by Eq. (11) for different $A$
parameters.

Figure 2: Initial time dependence of the survival probability $S(t)$ given by
Eq. (11) for parameter $A=0.3$ (broken line) in comparison with approximate
parabolic behaviour given by Eq. (12), with $\gamma=(1-A^{2})/(4\tau^{2})$
(solid line).

##  4 Some remarks

Usually, the memory kernel $K(t-s)$ can lead to non-Markovian dynamics
depending on the structure and time scale of the processes. Also a Markovian
stochastic processes (where the time correlation
$\langle\eta(t)\eta(s)\rangle\sim e^{-|t-s|/\tau}$ as is given by Eq. (6)) can
lead to non-Markovian dynamics of the system that is coupled to. In general,
Markovianity or non-Markovianity are not just features of the noise but of the
dynamics of the system coupled to the noise [17].

Survival probability given by Eq. (11) formally is the superposition

$S\left(t\right)=\omega_{a}\exp\left(-t/\tau_{a}\right)-\omega_{b}\exp\left(-t/\tau_{b}\right),$
(13)

with weights $\omega_{a}=\left(1+A\right)/\left(2A\right)$,
$\omega_{b}=\left(1-A\right)/\left(2A\right)$, and time constants
$\tau_{a}=\left(2\tau\right)/\left(1-A\right)$, and
$\tau_{b}=2\tau\left(1+A\right)$. For $t$ small enough we have

$S\left(t\right)\simeq\left(\omega_{a}-\omega_{b}\right)+\left(\frac{\omega_{b}}{\tau_{b}}-\frac{\omega_{a}}{\tau_{a}}\right)t+\left(\frac{\omega_{a}}{2\tau^{2}_{a}}-\frac{\omega_{b}}{2\tau^{2}_{b}}\right)t^{2}.$
(14)

From normalization, $S(t=0)=1$, we have $\omega_{a}-\omega_{b}=1$. For
$\omega_{a}/\tau_{a}=\omega_{b}/\tau_{b}$ and $\tau_{a}>\tau_{b}$ we observe
$t^{2}$-form od initial decay (as given by Eq. (12)). However, for
$\omega_{a}/\tau^{2}_{a}=\omega_{b}/\tau^{2}_{b}$ and $\tau_{a}<\tau_{b}$ we
have “normal” behaviour with $S(t)\simeq
1-(\omega_{b}/\tau_{b}-\omega_{a}/\tau_{a})t$.

The above discussion applies to $4\gamma\tau^{2}<1$ cases, nevertheless
$4\gamma\tau^{2}>1$ is not excluded. Differential equation (10) corresponds to
damped harmonic oscillator. For $4\gamma\tau^{2}>1$ the system oscillate with
amplitude gradually decreasing to zero (under-damped oscillator). The
frequency of oscillations is given by
$\omega^{2}=\gamma-\left(4\tau^{2}\right)^{-1}$ and amplitude decay with rate
$1/2\tau$. 555Exponentially damped oscillating modes in decay of unstable
quantum systems were widely discussed in [18, 19, 20].

##  5 Memory destroy independence

Consider the survival probability $S(t)\simeq 1-at^{2}$, as given by Eqs. (2)
and (12), where the time interval $[0,t]$ is interrupted by $n$ measurements
at times $t/n,2t/n,...,t$. The width of bins is $\delta=t/n$. Conditional
probability for $i-th$ bin is

$\displaystyle
S_{i}=S\left(t_{i}=i\delta|t_{i-1}=(i-1)\delta\right)=\frac{S(i\delta\cap(i-1)\delta)}{S((i-1)\delta)}=$
$\displaystyle=\frac{S(i\delta)}{S((i-1)\delta)}=\frac{1-a\delta^{2}i^{2}}{1-a\delta^{2}(i-1)^{2}}\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ (15)

and is not the same as $S(\delta)=1-a\delta^{2}$, contrary to the memory-less
exponential $S(t)$ for which $S_{i}=S(\delta)=e^{-\delta}$. For
$\delta\rightarrow 0$, the survival probability in each one of bins
$S_{i}\rightarrow 1$. Nevertheless, the survival probability for $n$
measurements

$S(t=n\delta)=\prod_{i=1}^{n}\frac{1-a\delta^{2}i^{2}}{1-a\delta^{2}(i-1)^{2}}=1-at^{2}$
(16)

and is not dependent on the number of measurements (in particular,
$S(n\delta)$ is the same when $n\rightarrow\infty$). Contrary to arguments
justifying the quantum Zeno effect, in classical system an unstable particle
observed continuously decay “normally”.

Exponential decay law, being memory-less, obeys for conditional survival
probability $S(t|t_{1})$ the relation $S(t=t_{1}+t_{2}|t_{1})=S(t_{2})$. Due
to independence we have $S(t_{1}+t_{2})=S(t_{1})S(t_{2})$ (in particular, such
assumption leads to Eq. (4) ). This is not valid for dependent variables, and

$S(t=t_{1}+t_{2})\neq S(t_{1})S(t_{2}).$ (17)

For dependent variables (due to noise correlation) we have chosen to use the
joint survival probability, in the form (in analogy to the power expansion of
multivariate normal distribution)

$S(t_{1},t_{2},...,t_{n})=1-\sum^{n}_{i=1}\sum^{n}_{j=1}t_{i}t_{j}C_{i,j},$
(18)

where
$C_{i,j}=\langle\eta(t_{i})\eta(t_{j})\rangle-\langle\eta(t_{i})\rangle\langle\eta(t_{j})\rangle$
and for noise $\eta(t)$ we can write $C_{i,j}=C$. Even for non-correlated
noise, $C_{i,j}=0$ for $i\neq j$, we have
$S(t_{1},t_{2},...,t_{n})=1-\sum^{n}_{i}t^{2}_{i}C_{i,i}\neq
S(\sum^{n}_{i}t_{t})$.

Interrupted the interval $[0,t]$ by $n$ measurements at times $\delta=t/n$,
the survival probability is

$S(\delta_{1},\delta_{2},...,\delta_{n})=1-C\sum^{n}_{i=1}\sum^{n}_{j=1}\left(\frac{t}{n}\right)^{2}=1-Ct^{2}$
(19)

and the Zeno paradox does not exist 666The absence of Zeno effect in dynamics
with classical evolution was considered recently in Ref. [21].

##  6 Summary

The dynamical evolution of statistical system is always influenced by their
environment which exhibits time-correlated random fluctuations which can lead
to non-Markovian dynamics [22, 23]. By studying time-correlated noise we have
shown how survival probability depends on the time scale of the noise
correlations. The time correlations in the noise field determine how fast the
survival probability converges to its exponential law behaviour.
Approximations, considering only one variable survival probability $S(t)$ do
not take into account all memory effects. The assumption of
$S\left(\sum_{i=1}^{n}t_{i}\right)=\prod_{i=1}^{n}S(t_{i})$, that leads to the
conclusion expressed by Eq. (4), is not valid in classical physics.
Independence of sequential measurements in quantum mechanics arises due to the
fact that, unlike in classical physics, measurements generally cannot be done
without changing the state of the measured system. In contrast, in classical
physics where measurements can be carried without affecting the state of the
system, there is no Zeno effect even for quadratic form of the survival
probabilities.

This research was supported by the National Science Centre, Poland (NCN) Grant
2020/39/O/ST2/00277 (MR). In preparation of this work we used the resources of
the Center for Computation and Computational Modeling of the Faculty of Exact
and Natural Sciences of the Jan Kochanowski University of Kielce.

##  Appendix A Formal solution of Eq. (5)

The formal solution of Eq. (5) can be obtained by means of Laplace transform,
namely:

$P\left(t\right)=\frac{1}{2\pi\imath}\int_{\sigma-\imath\infty}^{\sigma+\imath\infty}\frac{P_{0}+N\left(s\right)}{s+K\left(s\right)}e^{st}ds,$
(20)

where $P_{0}=P\left(t=0\right)$, and $K\left(s\right)$ and $N\left(s\right)$
denote the Laplace transforms of the memory kernel and the noise,
respectively. The integral is understood on a Bromwich contour, where $\sigma$
is a vertical contour in the complex plain chosen so that all singularities
are to the left of it.

We can evaluate the following averages:

$S\left(t\right)=\langle
P\left(t\right)\rangle=P_{0}\frac{1}{2\pi\imath}\int_{\sigma-\imath\infty}^{\sigma+\imath\infty}\frac{1}{s+K\left(s\right)}e^{st}ds,$
(21)

and

$\langle\left(P\left(t\right)-\langle
P\left(t\right)\rangle\right)^{2}\rangle=L^{-1}\left[\frac{\langle
N\left(s\right)N\left(s^{\prime}\right)\rangle}{\left(s+K\left(s\right)\right)\left(s^{\prime}+K\left(s^{\prime}\right)\right)}\right]\left(t,t^{\prime}\right),$
(22)

where $L^{-1}\left[h\right]\left(t,t^{\prime}\right)$ is the two-dimensional
inverse Laplece transform of $h\left(s,s^{\prime}\right)$ that depend on $t$
and $t^{\prime}$. The Laplace transform of the memory kernel, given by Eq. (6)
with noise correlator (8), is

$K\left(t\right)=\frac{\tau\gamma}{\tau s+1},$ (23)

and

$\langle
N\left(s\right)N\left(s^{\prime}\right)\rangle=\frac{\gamma\left(2+\tau s+\tau
s^{\prime}\right)}{\left(s+s^{\prime}\right)\left(1+\tau s\right)\left(1+\tau
s^{\prime}\right)}.$ (24)

From Eqs. (21) and (23), considering that the zeros of $s+K\left(s\right)$ are
the solutions of the equation $s+\tau s^{2}+\gamma=0$, we get Eq. (11) by a
straightforward application of the residues theorem.

## References

  * [1] H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford Univ. Press, Oxford, 2007).
  * [2] U. Weiss, Quantum Dissipative Systems, (World Scientific, Singapore, 2000).
  * [3] A. Rivas and S. F. Huelga, Open Quantum Systems. An Introduction, (Springer Briefs in Physics, Springer, 2011).
  * [4] L. Fonda, G. C. Ghirardi and A. Rimini, Rept. Prog. Phys. 41, 587-631 (1978) doi:10.1088/0034-4885/41/4/003
  * [5] J. Bogdanowicz, M. Pindor and R. Raczka, Found. Phys. 25, 833-849 (1995) doi:10.1007/BF02080566
  * [6] B. Misra and E. C. G. Sudarshan, J. Math. Phys. 18, 756 (1977) doi:10.1063/1.523304.
  * [7] K. Koshino and A. Shimizu, Phys. Rept. 412, 191-275 (2005) doi:10.1016/j.physrep.2005.03.001 [arXiv:quant-ph/0411145 [quant-ph]].
  * [8] P. Facchi and S. Pascazio, J. Phys. A 41, 493001 (2008) doi:10.1088/1751-8113/41/49/493001 [arXiv:0903.3297].
  * [9] P. Facchi, H. Nakazato, and S. Pascazio, Phys. Rev. Lett. 86, 2699 (2001) doi:10.1103/PhysRevLett.86.2699.
  * [10] C. Rothe, S. I. Hintschich and A. P. Monkman, Phys. Rev. Lett. 96, 163601 (2006) doi:10.1103/PhysRevLett.96.163601.
  * [11] G. Wilk and Z. Wlodarczyk, Phys. Lett. A 290, 55-58 (2001) doi:10.1016/S0375-9601(01)00645-4 [arXiv:hep-ph/0103114 [hep-ph]].
  * [12] W. M. Itano, Current Science 116, 201 (2019) doi: 10.18520/cs/v116/i2/201-204.
  * [13] D. Home and M. A. B. Whitaker, Annals of Physics 258, 237 (1997). doi: 10.1006/aphy.1997.5699.
  * [14] R. Zwanzig, Nonequilibrium Statistical Mechanics, (Oxford University Press, New York, 2001).
  * [15] M. Ruggieri, Pooja, J. Prakash and S. K. Das, Phys. Rev. D 106, no.3, 034032 (2022) doi:10.1103/PhysRevD.106.034032 [arXiv:2203.06712 [hep-ph]].
  * [16] P. Hanggi, K. E. Shuler, I. Oppennein, Physica A 107, 143 (1981) doi:10.1016/0378-4371(81)90028-5.
  * [17] M. M. Müller, S. Gherardini and F. Caruso, Sci. Rep. 6, 38650 (2016) doi:10.1038/srep38650.
  * [18] F. Giraldi, Eur. Phys. J. D 73, no.11, 239 (2019) doi:10.1140/epjd/e2019-100219-0 [arXiv:1903.01441 [quant-ph]].
  * [19] R. G. Winter, Phys. Rev. 123, 1503-1507 (1961) doi:10.1103/PhysRev.123.1503
  * [20] D. F. Ramírez Jiménez and N. G. Kelkar, Phys. Rev. A 104, no.2, 022214 (2021) doi:10.1103/PhysRevA.104.022214 [arXiv:2108.10957 [quant-ph]].
  * [21] J. López, L. Ares and A. Luis, J. Phys. A 53, no.37, 375306 (2020) doi:10.1088/1751-8121/abab53 [arXiv:2005.00265 [quant-ph]].
  * [22] C. Benedetti, M. Paris, and S. Maniscalco, Phys. Rev. A 89, 012114 (2014) doi:10.1103/PhysRevA.89.012114.
  * [23] C. Benedetti, F. Buscemi, P. Bordone, and M. Paris, Phys. Rev. A 93, 042313 (2016) doi:10.1103/PhysRevA.93.042313.