[2]Miguel A. Bessa

1]Department of Materials Science and Engineering, Delft University of
Technology, Mekelweg 2, Delft, 2628 CD, The Netherlands

2]School of Engineering, Brown University, 184 Hope St., Providence, RI 02912,
USA

# Continual learning for surface defect segmentation by subnetwork creation
and selection

Aleksandr Dekhovich<EMAIL_ADDRESS>[  [

###### Abstract

We introduce a new continual (or lifelong) learning algorithm called LDA-CP&S
that performs segmentation tasks without undergoing catastrophic forgetting.
The method is applied to two different surface defect segmentation problems
that are learned incrementally, i.e. providing data about one type of defect
at a time, while still being capable of predicting every defect that was seen
previously. Our method creates a defect-related subnetwork for each defect
type via iterative pruning and trains a classifier based on linear
discriminant analysis (LDA). At the inference stage, we first predict the
defect type with LDA and then predict the surface defects using the selected
subnetwork. We compare our method with other continual learning methods
showing a significant improvement – mean Intersection over Union better by a
factor of two when compared to existing methods on both datasets. Importantly,
our approach shows comparable results with joint training when all the
training data (all defects) are seen simultaneously111Code implementation is
available at: https://github.com/adekhovich/continual_defect_segmentation.

###### keywords:

Continual learning, Automatic vision inspection, Surface defect segmentation,
Linear discriminant analysis (LDA).

##  1 Introduction

Automatic defects inspection plays an important role in product quality
evaluation [1]. In the beginning of the field, the creation of meaningful
features to find defective regions was done manually [2, 3, 4, 5]. Although
classical machine learning methods have been proposed to identify images with
defective surfaces [6, 7, 8], recent advances in deep learning research have
led to an increase in performance [1]. Typically, there are three types of
tasks for defect inspection with neural networks – classification, detection
[9] and segmentation [10]. In the case of defect classification, transfer
learning helps to increase the network’s ability to detect defective surfaces
[11, 12]. For segmentation, most methods are based on the U-Net architecture
[13] taking advantage of convolutional layers that automatically extract
features from the images of the surfaces [14, 15, 16, 17]. Attention
mechanisms [18] employed in the model’s architecture can lead to even more
accurate predictions [19, 20].

The advent of deep learning models came with more data for training and
comparing these models in different real-life scenarios. For instance, after
[4] proposed their NEU-DET dataset with Hot Rolled Steel Strip Surface
defects, containing six types of defects, other groups collected datasets with
either different defect categories or a more significant number of defects,
e.g., GC10-DET [21] and X-SDD [22]. In segmentation literature, we can also
find examples of different categorizations of surface defects, e.g., the
Magnetic tile dataset [17] contains images of five types of defects together
with defect-free cases. As a final example, the dataset collected by [23] also
contains a large number of images but only has three types of defects.

Notwithstanding the increase in availability of datasets, there are many
instances where there are few types of defects in each dataset. This is a
natural occurrence in Engineering practice because many processes are not
amenable to high-throughput. Simultaneously, if new defects occur or if
another defect identification task with similar characteristics is
encountered, using the original dataset and neural network model while
considering new types of defects in similar (or even different) materials can
be invaluable. However, training the same neural network model on a new
dataset currently requires retraining it on all the data, even if the model
was already capable of detecting some types of defects. This happens because
deep learning models suffer from catastrophic forgetting [24, 25, 26]. In
conventional training, neural networks cannot learn new tasks without
forgetting old ones if the tasks are learned incrementally. Instead, the
continual learning field [27] aims to solve this type of problem where the
model receives data in batches (tasks) but aims to learn information
mitigating the forgetting issues.

We illustrate the impact of catastrophic forgetting on segmentation tasks in
Figure 1 by considering the defect segmentation dataset SD-saliency-900 [28].
This dataset consists of images with three types of defects: scratches,
patches and an inclusion. We illustrate this phenomenon by focusing on three
typical learning scenarios: 1) single-task training where each defect is
learned with a single network, meaning there are three networks in total; 2)
joint training where the model has access to the entire dataset at once; 3)
finetuning, in which the network learns to segment sequentially, adapting the
parameters for the new task, having them pretrained on previous ones.

Figure 1: Example of forgetting in the case of incremental learning of three
types of defects.

For all three learning scenarios, we quantify the segmentation performance via
the mean Intersection over Union (mIoU) score for every task after each
incremental step (Figure 1). We observe that finetuning on a new task leads to
a significant drop in performance for the previous task(s), as indicated by
the blue bars – a clear illustration that learning a sequence of tasks with a
single network leads to forgetting the previous tasks in the sequence
(catastrophic forgetting). However, forgetting does not occur in the case of
single- and joint-task training because the network is capable of learning
each of the defects separately without any pretraining, while also being
capable of learning all of them together. We also note that both single- and
joint-task training have comparable performance, despite a small decrease in
the latter case222As a short note, marginal improvements in performance
sometimes occur when changing the task order (investigated at the end of the
article). For example, the mIoU performance for the Scratches task improved by
0.13 points after learning the Inclusion task, but the improvement is small
compared to how much it degrades after learning the Patches task..

Therefore, we see that the ability to predict defects of previous types is
lost when training for a new type of defect, i.e. that is out-of-distribution.
The main objective of our work is to propose a continual learning algorithm
suitable for the surface segmentation problem. To the best of our knowledge,
this is the first work that develops a continual learning approach for surface
defect segmentation.

##  2 Continual learning

Overcoming the above-mentioned catastrophic forgetting requires deep learning
models to be trained in a continual (or lifelong) learning manner [29]. The
overwhelming majority of continual learning literature is dedicated to
classification tasks. In that context, three different categories have emerged
[27]: regularization-based [30, 31, 32], replay-based [33, 34, 9, 35] and
architectural-based methods [36, 37, 38]. In some cases, there are methods
that can fall into more than one category [39, 40]. Often these methods show
better performance but require more memory or have high computational cost
(e.g., extra memory buffer or architecture extension), which creates specific
challenges when deployed in real-life applications.

Regularization-based methods penalize parameters obtained on incremental step
$t-1$ from drastic changes while learning the task on incremental step $t$.
For example, SI [31], EWC [41] and MAS [32] employ total loss
$\mathcal{L}^{(t)}(x;\theta^{(t)})$ on incremental step $t$ that consists of
the loss computed for the current data, and a penalty term to prevent
forgetting:

$\mathcal{L}^{(t)}\big{(}x;\theta^{(t)}\big{)}=\mathcal{L}_{curr}\big{(}x;\theta^{(t)}\big{)}+\frac{\lambda}{2}\sum_{i=1}^{\\#params}\Omega_{i}\big{(}\theta_{i}^{(t)}-\theta_{i}^{(t-1)}\big{)}^{2},$
(1)

where $\mathcal{L}_{curr}\big{(}x;\theta^{(t)}\big{)}$ is a loss on the
current data,
$\sum_{i=1}^{\\#params}\Omega_{i}\big{(}\theta_{i}^{(t)}-\theta_{i}^{(t-1)}\big{)}^{2}$
is the penalty term, $\Omega_{i}$ is the cumulative importance for parameter
$i$, and $\theta^{(t-1)},\theta^{(t)}$ are network parameters at incremental
steps $t-1$ and $t$, respectively. Learning without forgetting (LwF) [30] aims
to mitigate forgetting by minimizing the cross-entropy between output
probabilities before and after the model is trained on a new task.

Replay-based (or rehearsal-based) approaches use a small fraction of data from
previous tasks and keep it in a fixed-size memory buffer. However, storing old
data in the buffer may not be allowed due to privacy issues [42]. Also, if the
model parameters were downloaded without the memory buffer, further model
training is not possible without forgetting. Therefore, in this work, we do
not focus on this type of methods.

Architectural approaches do manipulations with the network structure by
freezing and assigning some parameters to a specific task (e.g., PackNet [36])
or constantly growing the architecture increasing the expressivity of the
network (e.g., DEN [43]). However, if the model grows while learning a new
task, the final number of parameters is not bounded, leading to additional
computational costs. Alternatively, if the algorithm finds task-specific
parameters, e.g., via iterative pruning in CP&S [38] or pruning at
initialization in SupSup [44], the challenge lies on activating the correct
subnetwork during inference. This subnetwork selection in both CP&S and SupSup
requires a batch of test data of the task to be predicted, such that the
correct subnetwork (i.e. task ID) is identified. This may also be impractical
in real-life cases.

Literature on continual learning for semantic segmentation is scarcer. We can
find examples that adapt classification continual learning algorithms to
segmentation [45, 46], or some new approaches designed specifically for
segmentation [47, 48, 49]. Similar to the classification case, better results
are achieved by the methods that use a fixed-size memory buffer with samples
from old tasks to overcome forgetting [50, 51]. However, even though these
methods use old data (facilitating training), they still show significant
forgetting of the first tasks while performing well only on the last ones.

Continual learning also finds its application in industrial and manufacturing
cases. For example, MAS [32] was applied for product quality evaluation [52].
The approach clones the output head for previous tasks with the lowest loss on
the current data and uses this copy as initialization for a new task. The
weight transfer for the output layer, and MAS algorithm that penalizes
parameters from previous layers, show good performance for the considered
regression problem. Regularization-based methods have been examined for
anomaly detection in manufacturing process [53] and fault prediction in
lithium-ion batteries [54]. [55] developed an adaptive classification
framework based on continual learning to identify new unlabeled samples. The
proposed approach uses Mahalanobis distance and is employed to decide whether
a new batch of data belongs to the already seen defect type, or forms a new
one.

##  3 Proposed approach

We propose to take advantage of architectural methods that create task-
specific subnetworks for each task, eliminating the subnetwork selection
issue. As a base method, we consider Continual Prune-and-Select (CP&S) [38]
where we improve the subnetwork selection process by training a model for this
purpose, instead of having simple metric-based decision rules. In general, the
task-prediction problem is quite challenging in continual learning [56] and
can be seen as an out-of-distribution (OOD) detection problem [57]. The
difficulty arises from the presence of arbitrary classes in each task, leading
to cases where classes within each task may not be similar, while classes from
different tasks may have important similarities. This poses a challenge to
identify the task-ID and corresponding subnetwork, affecting the performance
of the continual learning model when the wrong subnetwork is selected.
Conversely, these methods have the advantage that when the correct subnetwork
is identified then there is no forgetting, which explains their state-of-the-
art performance in different image-classification datasets [38].

However, in contrast to image classification, every task in defect
segmentation problems consists of defects of only one type. This represents an
opportunity for architectural continual learning methods because we can train
a model that learns the distribution of each defect separately. To do so, we
use linear discriminant analysis trained on features extracted from a
pretrained convolutional neural network [58, 59]. For the segmentation model,
we use the U-Net architecture [13], in which we create task-specific
subnetworks via iterative pruning. As a pretrained feature extractor, we use
the EfficientNet-B5 architecture [60] pretrained on ImageNet-1000 [61].

Figure 2: An overview of the proposed method: (a) task ID (i.e., defect
type/subnetwork) prediction; (b) defects segmentation.

In Figure 2, we illustrate the inference stage of our approach, which consists
of two steps: (1) predicting the defect type (task ID) with LDA; and (2) using
a subnetwork that corresponds to the predicted defect to predict the
segmentation mask. Note that at the inference stage, defect type prediction
and defect mask prediction need to be done sequentially. Training for these
steps can be done in parallel and independently from each other. We call our
proposed approach LDA-CP&S since it uses the CP&S paradigm of creating
subnetworks during training, and it employs LDA for the subnetwork selection.

Referring again to Figure 1, we recall that the three separate models (single-
task grey bars) are capable of learning the defects slightly better than joint
training with all the tasks together (orange bars). This hints that having
task-specific parameters associated with only one task can even help the
learning process. At the same time, the shared parameters provide a transfer
learning effect between a new subnetwork and all the ones created before. Both
advantages can be exploited by LDA-CP&S.

Notwithstanding, our method could suffer a performance drop from two possible
sources: the pruning stage, and the LDA classification stage. The performance
reduction due to pruning may occur because some important parameters could be
deleted when creating additional space (free connections) for future tasks. In
addition, misclassification by LDA could result in signal routing through the
wrong subnetwork and consequently poor segmentation performance. In Section 4,
we show that these two sources of error are negligible compared to the
benefits of our approach. In the following subsections, we describe the
processes for subnetwork creation and LDA training.

###  3.1 Subnetwork creation

To create a subnetwork for the given task, we use NNrelief pruning algorithm
[62]. The approach evaluates the strength of the signal that propagates
through every connection/kernel. This pruning technique shows better sparsity
results than other connection/kernel-based pruning techniques [63, 64, 65].

For the set of $m_{l-1}$-channelled input samples
$\mathbf{X}^{l-1}=\\{\mathbf{x}^{l-1}_{1},\ldots,\mathbf{x}^{l-1}_{N}\\}$,
where
$\mathbf{x}^{(l-1)}_{k}=(x^{l-1}_{k1},\ldots,x^{l-1}_{km_{l-1}})\in\mathbb{R}^{m_{l-1}\times
h_{l-1}^{1}\times h^{2}_{l-1}}$ with $h^{1}_{l-1}$ and $h^{2}_{l-1}$ being the
height and width of feature maps for convolutional layer $l$. For every kernel
$\mathbf{K}^{l}_{1j},\mathbf{K}^{l}_{2j},\ldots,\mathbf{K}^{l}_{m_{l}j},\
\mathbf{K}^{l}_{ij}=(k^{l}_{ijqt})\in\mathbb{R}^{r_{l}\times r_{l}}$, $q\geq
1,\ r_{l}\geq t$, where $r_{l}$ is a kernel size, and for every bias
$b^{(l)}_{j}$ in filter $\mathbf{F}^{l}_{j}$, we define
$\mathbf{\hat{K}}^{l}_{ij}=\Big{(}\Big{\lvert}k^{l}_{ijqt}\Big{\rvert}\Big{)}$
as a matrix consisting of the absolute values of the matrix
$\mathbf{K}^{(l)}_{ij}$. Then we compute importance scores
$s^{l}_{ij},i\in\\{1,2,\ldots,m_{l}\\}$ of kernels $\mathbf{K}^{l}_{ij}$ as
follows:

$\displaystyle s^{l}_{ij}$
$\displaystyle=\frac{\frac{1}{N}\sum_{n=1}^{N}\Big{\lvert}\Big{\lvert}\mathbf{\hat{K}}^{l}_{ij}*\Big{\lvert}x^{l-1}_{ni}\Big{\rvert}\
\Big{\rvert}\Big{\rvert}_{F}}{S^{l}_{j}},$ (2)

where
$S^{l}_{j}=\sum_{i=1}^{m_{l-1}}\Big{(}\frac{1}{N}\sum_{n=1}^{N}\Big{\lvert}\Big{\lvert}\mathbf{\hat{K}}^{l}_{ij}*\Big{\lvert}x^{l-1}_{ni}\Big{\rvert}\
\Big{\rvert}\Big{\rvert}_{F}\Big{)}$ is the total importance score in filter
$\mathbf{F}^{l}_{j}$ of layer $l$, with $*$ indicating a convolution
operation, and where $\lvert\lvert\cdot\rvert\rvert_{F}$ is the Frobenius
norm.

The sketch of the algorithm for pruning filter $\mathbf{F}^{l}_{j}$ in a
convolutional layer $l$ can be described as follows:

  1. 1.

Choose $\alpha\in(0,1)$ – the amount of kernels’ importance that we want to
keep relative to the total importance of the kernels in the filter
$\mathbf{F}^{l}_{j}$.

  2. 2.

Compute importance scores $s_{ij}^{l}$ for all kernels in the filter
$\mathbf{F}^{l}_{j},\ i=1,\ldots,m_{l-1}$, using Eq. 2.

  3. 3.

Sort importance scores $s_{ij}^{l}$ for the filter $\mathbf{F}^{l}_{j}$.

  4. 4.

For the sorted importance scores $\hat{s}_{ij}^{l}$ find minimal $p\leq
m_{l-1}$ such that $\sum_{i=1}^{p}\hat{s}_{ij}^{l}\geq\alpha$.

  5. 5.

Prune kernels with the importance score $s_{ij}^{l}<\hat{s}_{pj}^{l}$ for all
$i\leq m_{l-1}$ and fixed $j$.

Figure 3: An overview of Continual Prune-and-Select (CP&S): an example with
two tasks.

Overall, NNrelief finds kernels that propagate on average the lowest signal
according to the Frobenious norm and prune these kernels. As the outcome of
the procedure, we obtain a subnetwork (sub-U-Net) that predicts the defect for
only one type of defects. Then we fix all parameters that are assigned to this
subnetwork and do not update them anymore. When the network receives a new
task with a new type of defect, CP&S finds a subnetwork for this task within
the main U-Net, using the parameters assigned to the previous tasks, but
without updating them. Algorithm 1 illustrates the pseudocode for CP&S and
Figure 3 illustrated the method:

1:network $\mathcal{N}$, dataset $\\{\textbf{X}^{t}\\}_{t=1}^{T}$. Initialize
learning parameters (learning rate, weight decay, number of epochs, etc. ),
pruning parameters (for NNrelief algorithm: $\alpha$ and $num\\_iters$).

2:for $t=1,2,\ldots,T$ do

3: $\mathcal{N}^{(t)}\leftarrow\mathcal{N}$

4: for $it=1,2,\ldots,num\\_iters$ do $\triangleright$ repeat pruning

5:
$\mathcal{N}^{(t)}\leftarrow\text{Pruning}(\mathcal{N}^{(t)},\mathbf{X}^{(t)},\alpha)$
$\triangleright$ pruning step: NNrelief

6: Retrain subnetworks $\mathcal{N}^{(t)}$ $\triangleright$ retraining step

7: end for

8: Freeze parameters $w\in\mathcal{N}^{(t)}$ and never update them

9:end for

10:network $\mathcal{N}$ that learned tasks $1,2,\ldots,T$.

Algorithm 1 Pseudocode for CP&S training procedure

###  3.2 Subnetwork prediction (or selection)

To predict the task ID (type of defect) at the inference stage, we propose to
use linear discriminant analysis (LDA). In this subsection, we describe the
training procedure for LDA. In LDA, it is assumed that all classes have class
means $\mathbf{\mu}^{(1)},\mathbf{\mu}^{(2)},\ldots,\mathbf{\mu}^{(T)}$ and
share the same covariance matrix $\mathbf{\Sigma}$. However, in continual
learning, we do not have access to all tasks at the same time, but only task
$t$. Therefore, the covariance matrix needs to be updated online with respect
to the new data batch.

Let us denote a new given task as
$\mathbf{X}^{(t)}=\\{x_{1}^{(t)},x_{2}^{(t)},\dots,x_{N_{t}}^{(t)}\\}$.
Following streaming LDA (SLDA) strategy [59], we use a feature extractor
$\mathcal{F}$ pretrained on ImageNet-1000 to obtain low-dimensional data
representation
$\mathbf{Z}^{(t)}:=\\{z^{(t)}_{1},z^{(t)}_{2},\ldots,z^{(t)}_{N_{t}}\\}$,
$z^{(t)}_{i}=\mathcal{F}(x_{i}^{(t)})\in\mathbb{R}^{d}$. Then we can compute
the class mean $\mathbf{\mu}^{(t)}\in\mathbb{R}^{d}$ and update the shared
covariance matrix $\mathbf{\Sigma}^{(1:t)}\in\mathbb{R}^{d\times d}$ after
incremental step $t$ as follows [66]:

$\displaystyle\mathbf{\mu}^{(t)}$
$\displaystyle=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}z_{i}^{(t)}$ (3)
$\displaystyle\mathbf{\Sigma}^{(1:t)}$
$\displaystyle=\frac{(t-1)\mathbf{\Sigma}^{(1:t-1)}+\mathbf{\Delta}^{(t)}}{t},$
(4)

where
$\mathbf{\Delta}^{(t)}=\frac{(t-1)(Z^{(t)}-\mu^{(t)})(Z^{(t)}-\mu^{(t)})^{\mathsf{T}}}{t}$
and
$(Z^{(t)}-\mu^{(t)}):=(z_{1}^{(t)}-\mu^{(t)},z_{2}^{(t)}-\mu^{(t)},\ldots,z_{N_{t}}^{(t)}-\mu^{(t)})\in\mathbb{R}^{d\times
N_{t}}$. In SLDA, the regularized version of LDA is implemented by applying
shrinkage regularization to covariance matrix:
$\mathbf{\Lambda}^{(1:t)}=[(1-\varepsilon)\mathbf{\Sigma}^{(1:t)}+\varepsilon\mathbf{I}]^{-1}$,
where $\mathbf{I}$ is an identity matrix of the corresponding dimension.

At the inference stage, after learning all class means $\mathbf{\mu}^{(t)},\
t=1,2,\ldots,T$, and shared covariance matrix $\mathbf{\Sigma}^{(1:T)}$ (and
$\mathbf{\Lambda}^{(1:t)}$ as a result), we can make a prediction for a new
test sample $x$ as follows:

$c=\underset{i=1,2,\ldots,T}{\text{argmax
}}(\mathbf{W}\mathcal{F}(x)+\mathbf{b})_{i},$ (5)

where $\mathbf{W}=\mathbf{M}^{(1:T)}\mathbf{\Lambda}^{(1:T)}$, rows of
$\mathbf{M}^{(1:T)}$ are mean vectors $\mathbf{\mu}^{(t)}\ (t=1,2,\ldots,T)$,
and $\mathbf{b}_{i}=-\frac{1}{2}\mu^{(i)}\mathbf{\Lambda}^{(1:T)}\mu^{(i)}$.

Unlike previous task prediction strategies [44, 67], with LDA we can predict
the task ID with a single test sample, rather than with a batch of samples,
representing an important advantage. This is possible because each task
consists of defects of the same type and can be described well by a normal
distribution with class means $\mathbf{\mu}^{(t)}$ and common covariance
matrix $\mathbf{\Sigma}^{(1:T)}$.

##  4 Experiments and results

We evaluate our LDA-CP&S approach on the SD-saliency-900 [28] and Magnetic
tile defects [17] datasets, comparing with the following scenarios:

  * •

joint training: the model has access to all data at each incremental step.
This case is an upper bound for rehearsal-based methods.

  * •

finetuning: the model is trained at each incremental step $t$ without
preventing forgetting, i.e., we finetune the model to a new task $t$ that is
pretrained on previous tasks $1,2,\ldots,t-1$, inevitably causing forgetting
of previous tasks because the network parameters (weights and biases) are
updated for task $t$.

  * •

Regularization-based continual learning methods: LwF, MAS that penalize
important parameters from changing (see Section 2, Eq. 1), in an attempt to
alleviate forgetting.

We do not consider rehearsal-based approaches that replay a small portion of
data from previous tasks while learning a new one because our premise is that
old data is not available and should not be used. Furthermore, our comparative
investigation of the proposed LDA-CP&S method with others includes the joint
training strategy, which is an upper bound for rehearsal-based methods, where
all data is available at each incremental step. Therefore, if we show that
LDA-CP&S performs similarly to joint training, there is no need to consider
rehearsal-based continual learning methods.

As performance metrics, we follow other segmentation works and use the mean
Intersection over Union score. For ground truth $Y$ and prediction $\hat{Y}$,
IoU score is computed as follows:

$\displaystyle\text{IoU}(\hat{Y},Y)=\frac{\lvert Y\cap\hat{Y}\rvert}{\lvert
Y\cup\hat{Y}\rvert}.$ (6)

To train the model, we use IoU loss which leads to better performance in our
experiments than other losses, e.g., Tversky loss [68] and Focal loss [69].
However, it is worth noting that the difference in IoU scores between models
trained with different loss functions is not significant. The IoU loss is
computed as follows:

$\displaystyle\text{IoUloss}(\hat{P},Y)=1-\frac{\sum_{i=1}^{H}\sum_{j=1}^{W}\hat{P}_{ij}\cdot
Y_{ij}+\varepsilon}{\sum_{i=1}^{H}\sum_{j=1}^{W}\hat{P}_{ij}+Y_{ij}-\hat{P}_{ij}\cdot
Y_{ij}+\varepsilon},$ (7)

where $\hat{P}_{ij}\in[0,1]^{H\times W}$ are the output probabilities, and $H$
and $W$ are the height and width of the output image, $\varepsilon$ is a
smoothing parameter.

###  4.1 SD-saliency-900 dataset

In the case of the SD-saliency-900 dataset, we consider a smaller version of
U-Net with 16, 32, 64 and 128 in the encoder block and 256 channels in the
bottleneck because it consists of only three types of defects – Scratches,
Patches and Inclusion – with 300 images per defect. The original size of the
images is $200\times 200$ but we resize the images to $224\times 224$ to make
them acceptable for U-Net. We train the segmentation model for 70 epochs with
8 images in a batch, using Adam [70] optimizer and learning rate $0.001$.
During the pruning stage, we use $\alpha=0.9$ and $3$ pruning iterations. More
details about the influence of hyperparameters on the results are shown in
Section 4.3. As it is common in continual learning literature [71], we
consider different task orderings in our experiments. We can construct six
task orderings for the current dataset (e.g., Patches $\to$ Scratches $\to$
Inclusion).

Table 1: Classification accuracy (%) for SD-saliency-900 dataset. The numbers are averaged over all six orderings. | Scratches | Patches | Inclusion |  | Average
---|---|---|---|---|---
accuracy (%) | 98.33 | 100 | 100 |  | 99.44
|  |  |  |  |

First, we have to make sure that LDA can accurately predict the defect type in
an incremental manner. Table 1 illustrates the classifier’s accuracy for each
defect averaged over all six orders. We can observe the high performance of
LDA, misclassifying only a few images from the Scratches dataset. Since 60
images were selected to test each defect type, the prediction error presented
corresponds to only 1 misclassified image.

Figure 4: IoU score after every incremental step for SD-saliency-900 dataset.
The results are presented for all six possible defect orderings.

In Figure 4, we show mIoU score after every incremental step for every task
order. Regularization-based methods only slightly outperform finetuning
strategy, while our LDA-CP&S shows comparable results to joint training. Poor
performance of the regularization methods can be explained by the lack of a
task-specific output layer, which is present in classification network
architectures as a classification head. Therefore MAS and LwF update all
parameters but change them slightly less than finetuning. On the contrary,
LDA-CP&S creates fixed task-specific subnetworks that can overlap and transfer
knowledge between each other. Since LDA predicts the defect type (i.e.,
subnetwork) well at the inference stage, we almost do not have any losses in
segmentation performance. We also do not observe network saturation, i.e., the
situation when the model does not have enough free space to learn a new task,
even though we use a smaller version of U-Net.

###  4.2 Magnetic tile defects dataset

Magnetic tile defects dataset [17] contains five types of defects, namely
Blowhole, Break, Crack, Fray and Uneven, and images that are free from defects
(Free). In this work, we consider only images with defects, i.e., five
classes. Since the number of defects is higher in this case, we use a U-Net of
the original size with 64, 128, 256, and 512 in the encoder block and 1024
channels in the bottleneck. All images in the dataset have different image
sizes and, therefore, we resize them to $224\times 224$. For every defect, we
randomly select 80% images for training and the rest for testing. U-Net is
trained for 150 epochs with 8 images in a batch, using Adam optimizer and
learning rate $0.0001$. Since the total number of possible task orderings is
quite large ($5!=120$), we consider only five of them at random and we do not
have reason to believe that the final performance would be very different when
choosing other orderings:

  * •

Blowhole $\to$ Break $\to$ Crack $\to$ Fray $\to$ Uneven;

  * •

Break $\to$ Uneven $\to$ Fray $\to$ Crack $\to$ Blowhole;

  * •

Crack $\to$ Blowhole $\to$ Break $\to$ Uneven $\to$ Fray;

  * •

Fray $\to$ Crack $\to$ Uneven $\to$ Blowhole $\to$ Break;

  * •

Uneven $\to$ Fray $\to$ Blowhole $\to$ Break $\to$ Crack,

where each defect type appears exactly once for each ordering.

One of the main difficulties with this dataset is class imbalance. Table 2
presents LDA accuracy with the same feature extractor considered in the
previous example: the pretrained EfficientNet-B5 architecture. Overall, our
classification model is able to identify correctly four out of five types of
defects, having some difficulties with the Fray sub-dataset that contains the
smallest number of images. The only mistake was done in the Fray sub-dataset
where we have only 7 test images, meaning that only one image is classified
wrongly.

Table 2: Classification accuracy (%) for Magnetic tile dataset and the total size of the dataset. The numbers for accuracy are averaged over all five orderings. | Blowhole | Break | Crack | Fray | Uneven |  | Average
---|---|---|---|---|---|---|---
# train (test) images | 92 (23) | 92 (22) | 68 (17) | 25 (7) | 72 (21) |  | N/A
test accuracy (%) | 100 | 100 | 100 | 85.71 | 100 |  | 98.75
|  |  |  |  |  |  |

We would like to highlight the necessity of pretraining the network that
extracts features for LDA. The pretrained EfficientNet-B5 produces lower
dimensional embeddings that can be used for training and classification with
an accuracy of 98.75%, misclassifying only one test image. Meanwhile, if we
were to consider a feature extractor with random parameters it would compress
the input images in such a way that the LDA classifier would only achieve
16.24% of accuracy.

Figure 5: IoU score after every incremental step for Magnetic tile datasets.
The results are shown for all five selected defect orderings.

In Figure 5, we present the mIoU score after every incremental step, comparing
our LDA-CP&S with other continual learning methods. As we saw in the previous
example, regularization-based methods do not handle this type of segmentation
problem well. The tasks that we constructed from the Magnetic tile dataset can
be quite dissimilar having significant differences in defects areas.
Therefore, by updating all the parameters without having task-specific ones,
regularization-based approaches are only slightly better than simple
finetuning where no anti-forgetting measures are considered. In contrast, our
LDA-CP&S creates task-specific parameters for each defect, fixing the values
of the parameters once they are assigned to a subnetwork (i.e., defect type or
task ID). This allows LDA-CP&S to deal with sequences of tasks as well as
joint training, which is very encouraging because joint training is a
performance upper bound since all the data is available at each incremental
step.

We also investigated how the mIoU score changes for every task after each
incremental step. In Figure 6, we consider one of the task orderings: Fray
$\to$ Crack $\to$ Uneven $\to$ Blowhole $\to$ Break. The figure clearly shows
the advantage of our algorithm over regularization-based ones because they are
heavily dependent on the similarity of the tasks in the order. For example,
learning the Break sub-dataset (the last incremental step) improves
performance on Fray and Crack sub-datasets compared to the previous
incremental step for MAS, LwF and finetuning strategies. However, Uneven is
totally forgotten after the network is trained on the Blowhole sub-dataset.

Figure 6: IoU score after every incremental step for one of the defects
orderings from the Magnetic tile dataset.

On the contrary, LDA-CP&S does not forget previous tasks and is still able to
learn new ones even having fewer free parameters. It has comparable
performance with a single-task scenario, where a separate U-Net is trained for
every task. Also, we observe that task-wise performance is almost the same as
for joint training, meaning the subnetwork overlaps provide enough knowledge
transfer to learn a new task.

Figure 7 illustrates the model output in every learning scenario. We observe
that regularization-based methods and finetuning cannot capture the defects of
the first tasks in the sequence, while our LDA-CP&S finds defects’ segments
close to the joint training.

Figure 7: Visualization of defects prediction for the considered scenarios on
Fray $\to$ Crack $\to$ Uneven $\to$ Blowhole $\to$ Break task ordering.

###  4.3 Hyperparameters choice

The choice of hyperparameters for pruning has a significant impact on
subnetwork sparsity and, as a result, performance. In this subsection, we
compare different options for the pruning hyperparameter $\alpha$ and the
number of pruning iterations. A lower number of $\alpha$ and a higher number
of pruning iterations lead to higher sparsity (more free connections to learn
future tasks) but may cause lower segmentation performance. Also, the values
for hyperparameters depend on the length of task sequences. In our work, we
pre-define these hyperparameters at the beginning and do not change them
during the training process.

Figure 8: Comparison between different pairs of hyperparameters for the
pruning step in our LDA-CP&S method on SD-saliency-900 dataset (left) and
Magnetic tile dataset (right). The results after each incremental step
averaged over the number of considered task orderings.

Figure 8 illustrates how different pairs of hyperparameters affect the
training process for our approach. For both datasets, we clearly see that the
network starts to saturate if pruning is not aggressive enough (e.g.,
$\alpha=0.95$ where most of the signal is conserved) because the network does
not have enough free parameters for new tasks. In the case of the SD-
saliency-900 dataset, we can also observe the trade-off between sparsity and
mIoU score: with $\alpha=0.9$ it is clear that pruning the network twice leads
to better performance than doing it three times, as the subnetwork that
results is less expressive (has fewer parameters). The results on the Magnetic
tile dataset show the trade-off between learning the first tasks and the last
ones: if we prune the network twice, $\alpha=0.9$ leads to better performance
if there are no more than three tasks, while $\alpha=0.85$ is better suitable
for longer task sequences.

##  5 Conclusion

We believe smart monitoring systems should quickly adapt to new tasks without
a dramatic drop in performance on previously learned ones. However, this is
not the case based on the current state of the literature on surface defect
inspection. Thus, there is a need for continual learning of deep neural
networks for automatic surface defect segmentation such that product quality
assessment is improved. By training deep learning models incrementally, we
show that we can accumulate all the learned information without retraining
when a new task comes to the network. In addition, we do not need to store
data for retraining, which can be either not allowed or not possible due to
the (lack of) availability of old datasets.

The LDA-CP&S method that we propose successfully learns to segment the defects
incrementally, without any forgetting, using only the data that is given at
the current time step. Meanwhile, other methods that do not use data from
previous tasks fail to remember all tasks, exhibiting considerable forgetting
in segmenting previously seen defects. Overall, the performance of LDA-CP&S is
more than two times higher in terms of mean Intersection over Union score for
the two datasets considered herein when compared to other continual learning
methods. Moreover, it is comparable with joint training where the model has
access to all the data observed up to the current incremental step.

## References

  *   * Prunella et al. [2023] Prunella, M., Scardigno, R.M., Buongiorno, D., Brunetti, A., Longo, N., Carli, R., Dotoli, M., Bevilacqua, V.: Deep learning for automatic vision-based recognition of industrial surface defects: a survey. IEEE Access (2023)
  * Ojala et al. [2002] Ojala, T., Pietikainen, M., Maenpaa, T.: Multiresolution gray-scale and rotation invariant texture classification with local binary patterns. IEEE Transactions on pattern analysis and machine intelligence 24(7), 971–987 (2002)
  * Chao and Tsai [2008] Chao, S.-M., Tsai, D.-M.: An anisotropic diffusion-based defect detection for low-contrast glass substrates. Image and Vision Computing 26(2), 187–200 (2008)
  * Song and Yan [2013] Song, K., Yan, Y.: A noise robust method based on completed local binary patterns for hot-rolled steel strip surface defects. Applied Surface Science 285, 858–864 (2013)
  * Jeon et al. [2014] Jeon, Y.-J., Choi, D.-c., Lee, S.J., Yun, J.P., Kim, S.W.: Defect detection for corner cracks in steel billets using a wavelet reconstruction method. JOSA A 31(2), 227–237 (2014)
  * Jia et al. [2004] Jia, H., Murphey, Y.L., Shi, J., Chang, T.-S.: An intelligent real-time vision system for surface defect detection. In: Proceedings of the 17th International Conference on Pattern Recognition, 2004. ICPR 2004., vol. 3, pp. 239–242 (2004). IEEE
  * Agarwal et al. [2011] Agarwal, K., Shivpuri, R., Zhu, Y., Chang, T.-S., Huang, H.: Process knowledge based multi-class support vector classification (pk-msvm) approach for surface defects in hot rolling. Expert Systems with Applications 38(6), 7251–7262 (2011)
  * Shanmugamani et al. [2015] Shanmugamani, R., Sadique, M., Ramamoorthy, B.: Detection and classification of surface defects of gun barrels using computer vision and machine learning. Measurement 60, 222–230 (2015)
  * He et al. [2019] He, Y., Song, K., Meng, Q., Yan, Y.: An end-to-end steel surface defect detection approach via fusing multiple hierarchical features. IEEE transactions on instrumentation and measurement 69(4), 1493–1504 (2019)
  * Tabernik et al. [2020] Tabernik, D., Šela, S., Skvarč, J., Skočaj, D.: Segmentation-based deep-learning approach for surface-defect detection. Journal of Intelligent Manufacturing 31(3), 759–776 (2020)
  * Aslam et al. [2020] Aslam, M., Khan, T.M., Naqvi, S.S., Holmes, G., Naffa, R.: Ensemble convolutional neural networks with knowledge transfer for leather defect classification in industrial settings. IEEE Access 8, 198600–198614 (2020)
  * Wu and Lv [2021] Wu, H., Lv, Q.: Hot-rolled steel strip surface inspection based on transfer learning model. Journal of Sensors 2021, 1–8 (2021)
  * Ronneberger et al. [2015] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, pp. 234–241 (2015). Springer
  * He et al. [2019] He, D., Xu, K., Zhou, P.: Defect detection of hot rolled steels with a new object detection framework called classification priority network. Computers & Industrial Engineering 128, 290–297 (2019)
  * Song et al. [2020] Song, G., Song, K., Yan, Y.: Edrnet: Encoder–decoder residual network for salient object detection of strip steel surface defects. IEEE Transactions on Instrumentation and Measurement 69(12), 9709–9719 (2020)
  * Hao et al. [2021] Hao, R., Lu, B., Cheng, Y., Li, X., Huang, B.: A steel surface defect inspection approach towards smart industrial monitoring. Journal of Intelligent Manufacturing 32, 1833–1843 (2021)
  * Huang et al. [2020] Huang, Y., Qiu, C., Yuan, K.: Surface defect saliency of magnetic tile. The Visual Computer 36, 85–96 (2020)
  * Vaswani et al. [2017] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, Ł., Polosukhin, I.: Attention is all you need. Advances in neural information processing systems 30 (2017)
  * Pan and Zhang [2022] Pan, Y., Zhang, L.: Dual attention deep learning network for automatic steel surface defect segmentation. Computer-Aided Civil and Infrastructure Engineering 37(11), 1468–1487 (2022)
  * Üzen et al. [2022] Üzen, H., Türkoğlu, M., Yanikoglu, B., Hanbay, D.: Swin-mfinet: Swin transformer based multi-feature integration network for detection of pixel-level surface defects. Expert Systems with Applications 209, 118269 (2022)
  * Lv et al. [2020] Lv, X., Duan, F., Jiang, J.-j., Fu, X., Gan, L.: Deep metallic surface defect detection: The new benchmark and detection network. Sensors 20(6), 1562 (2020)
  * Feng et al. [2021] Feng, X., Gao, X., Luo, L.: X-sdd: A new benchmark for hot rolled steel strip surface defects detection. Symmetry 13(4), 706 (2021)
  * Liu and Ye [2022] Liu, T., Ye, W.: A semi-supervised learning method for surface defect classification of magnetic tiles. Machine Vision and Applications 33(2), 35 (2022)
  * French [1999] French, R.M.: Catastrophic forgetting in connectionist networks. Trends in cognitive sciences 3(4), 128–135 (1999)
  * Goodfellow et al. [2013] Goodfellow, I.J., Mirza, M., Xiao, D., Courville, A., Bengio, Y.: An empirical investigation of catastrophic forgetting in gradient-based neural networks. arXiv preprint arXiv:1312.6211 (2013)
  * Coop et al. [2013] Coop, R., Mishtal, A., Arel, I.: Ensemble learning in fixed expansion layer networks for mitigating catastrophic forgetting. IEEE transactions on neural networks and learning systems 24(10), 1623–1634 (2013)
  * De Lange et al. [2021] De Lange, M., Aljundi, R., Masana, M., Parisot, S., Jia, X., Leonardis, A., Slabaugh, G., Tuytelaars, T.: A continual learning survey: Defying forgetting in classification tasks. IEEE transactions on pattern analysis and machine intelligence 44(7), 3366–3385 (2021)
  * Song et al. [2020] Song, G., Song, K., Yan, Y.: Saliency detection for strip steel surface defects using multiple constraints and improved texture features. Optics and Lasers in Engineering 128, 106000 (2020)
  * Thrun [1998] Thrun, S.: Lifelong learning algorithms. In: Learning to Learn, pp. 181–209. Springer, ??? (1998)
  * Li and Hoiem [2017] Li, Z., Hoiem, D.: Learning without forgetting. IEEE transactions on pattern analysis and machine intelligence 40(12), 2935–2947 (2017)
  * Zenke et al. [2017] Zenke, F., Poole, B., Ganguli, S.: Continual learning through synaptic intelligence. In: International Conference on Machine Learning, pp. 3987–3995 (2017). PMLR
  * Aljundi et al. [2018] Aljundi, R., Babiloni, F., Elhoseiny, M., Rohrbach, M., Tuytelaars, T.: Memory aware synapses: Learning what (not) to forget. In: Proceedings of the European Conference on Computer Vision (ECCV), pp. 139–154 (2018)
  * Rebuffi et al. [2017] Rebuffi, S.-A., Kolesnikov, A., Sperl, G., Lampert, C.H.: icarl: Incremental classifier and representation learning. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 2001–2010 (2017)
  * Castro et al. [2018] Castro, F.M., Marín-Jiménez, M.J., Guil, N., Schmid, C., Alahari, K.: End-to-end incremental learning. In: Proceedings of the European Conference on Computer Vision (ECCV), pp. 233–248 (2018)
  * Douillard et al. [2020] Douillard, A., Cord, M., Ollion, C., Robert, T., Valle, E.: Podnet: Pooled outputs distillation for small-tasks incremental learning. In: Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XX 16, pp. 86–102 (2020). Springer
  * Mallya and Lazebnik [2018] Mallya, A., Lazebnik, S.: Packnet: Adding multiple tasks to a single network by iterative pruning. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 7765–7773 (2018)
  * Sokar et al. [2022] Sokar, G., Mocanu, D.C., Pechenizkiy, M.: Avoiding forgetting and allowing forward transfer in continual learning via sparse networks. In: Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pp. 85–101 (2022). Springer
  * Dekhovich et al. [2023] Dekhovich, A., Tax, D.M., Sluiter, M.H., Bessa, M.A.: Continual prune-and-select: class-incremental learning with specialized subnetworks. Applied Intelligence, 1–16 (2023)
  * Yan et al. [2021] Yan, S., Xie, J., He, X.: Der: Dynamically expandable representation for class incremental learning. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 3014–3023 (2021)
  * Wang et al. [2022] Wang, F.-Y., Zhou, D.-W., Ye, H.-J., Zhan, D.-C.: Foster: Feature boosting and compression for class-incremental learning. In: European Conference on Computer Vision, pp. 398–414 (2022). Springer
  * Kirkpatrick et al. [2017] Kirkpatrick, J., Pascanu, R., Rabinowitz, N., Veness, J., Desjardins, G., Rusu, A.A., Milan, K., Quan, J., Ramalho, T., Grabska-Barwinska, A., et al.: Overcoming catastrophic forgetting in neural networks. Proceedings of the national academy of sciences 114(13), 3521–3526 (2017)
  * Zhang et al. [2020] Zhang, J., Zhang, J., Ghosh, S., Li, D., Tasci, S., Heck, L., Zhang, H., Kuo, C.-C.J.: Class-incremental learning via deep model consolidation. In: Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pp. 1131–1140 (2020)
  * Yoon et al. [2017] Yoon, J., Yang, E., Lee, J., Hwang, S.J.: Lifelong learning with dynamically expandable networks. arXiv preprint arXiv:1708.01547 (2017)
  * Wortsman et al. [2020] Wortsman, M., Ramanujan, V., Liu, R., Kembhavi, A., Rastegari, M., Yosinski, J., Farhadi, A.: Supermasks in superposition. Advances in Neural Information Processing Systems 33, 15173–15184 (2020)
  * Baweja et al. [2018] Baweja, C., Glocker, B., Kamnitsas, K.: Towards continual learning in medical imaging. In: Medical Imaging Meets NIPS Workshop, 32nd Conference on Neural Information Processing Systems (NIPS) (2018)
  * van Garderen et al. [2019] Garderen, K., Voort, S., Incekara, F., Smits, M., Klein, S.: Towards continuous learning for glioma segmentation with elastic weight consolidation. In: International Conference on Medical Imaging with Deep Learning –Extended Abstract Track (2019)
  * Klingner et al. [2020] Klingner, M., Bär, A., Donn, P., Fingscheidt, T.: Class-incremental learning for semantic segmentation re-using neither old data nor old labels. In: 2020 IEEE 23rd International Conference on Intelligent Transportation Systems (ITSC), pp. 1–8 (2020). IEEE
  * Douillard et al. [2021] Douillard, A., Chen, Y., Dapogny, A., Cord, M.: Plop: Learning without forgetting for continual semantic segmentation. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 4040–4050 (2021)
  * Yan et al. [2021] Yan, S., Zhou, J., Xie, J., Zhang, S., He, X.: An em framework for online incremental learning of semantic segmentation. In: Proceedings of the 29th ACM International Conference on Multimedia, pp. 3052–3060 (2021)
  * Cha et al. [2021] Cha, S., Yoo, Y., Moon, T., et al.: Ssul: Semantic segmentation with unknown label for exemplar-based class-incremental learning. Advances in neural information processing systems 34, 10919–10930 (2021)
  * Qiu et al. [2023] Qiu, Y., Shen, Y., Sun, Z., Zheng, Y., Chang, X., Zheng, W., Wang, R.: Sats: Self-attention transfer for continual semantic segmentation. Pattern Recognition 138, 109383 (2023)
  * Tercan et al. [2022] Tercan, H., Deibert, P., Meisen, T.: Continual learning of neural networks for quality prediction in production using memory aware synapses and weight transfer. Journal of Intelligent Manufacturing 33(1), 283–292 (2022)
  * Maschler et al. [2021] Maschler, B., Pham, T.T.H., Weyrich, M.: Regularization-based continual learning for anomaly detection in discrete manufacturing. Procedia CIRP 104, 452–457 (2021)
  * Maschler et al. [2022] Maschler, B., Tatiyosyan, S., Weyrich, M.: Regularization-based continual learning for fault prediction in lithium-ion batteries. Procedia CIRP 112, 513–518 (2022)
  * Sun et al. [2022] Sun, W., Kontar, R.A., Jin, J., Chang, T.-S.: A continual learning framework for adaptive defect classification and inspection. arXiv preprint arXiv:2203.08796 (2022)
  * Kim et al. [2020] Kim, E.S., Kim, J.U., Lee, S., Moon, S.-K., Ro, Y.M.: Class incremental learning with task-selection. In: 2020 IEEE International Conference on Image Processing (ICIP), pp. 1846–1850 (2020). IEEE
  * Kim et al. [2022] Kim, G., Xiao, C., Konishi, T., Ke, Z., Liu, B.: A theoretical study on solving continual learning. Advances in Neural Information Processing Systems 35, 5065–5079 (2022)
  * Dorfer et al. [2015] Dorfer, M., Kelz, R., Widmer, G.: Deep linear discriminant analysis. arXiv preprint arXiv:1511.04707 (2015)
  * Hayes and Kanan [2020] Hayes, T.L., Kanan, C.: Lifelong machine learning with deep streaming linear discriminant analysis. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, pp. 220–221 (2020)
  * Tan and Le [2019] Tan, M., Le, Q.: Efficientnet: Rethinking model scaling for convolutional neural networks. In: International Conference on Machine Learning, pp. 6105–6114 (2019). PMLR
  * Deng et al. [2009] Deng, J., Dong, W., Socher, R., Li, L.-J., Li, K., Fei-Fei, L.: Imagenet: A large-scale hierarchical image database. In: 2009 IEEE Conference on Computer Vision and Pattern Recognition, pp. 248–255 (2009). Ieee
  * Dekhovich et al. [2021] Dekhovich, A., Tax, D.M., Sluiter, M.H., Bessa, M.A.: Neural network relief: a pruning algorithm based on neural activity. arXiv preprint arXiv:2109.10795 (2021)
  * Han et al. [2015] Han, S., Pool, J., Tran, J., Dally, W.: Learning both weights and connections for efficient neural network. Advances in neural information processing systems 28 (2015)
  * Li et al. [2016] Li, H., Kadav, A., Durdanovic, I., Samet, H., Graf, H.P.: Pruning filters for efficient convnets. arXiv preprint arXiv:1608.08710 (2016)
  * Lee et al. [2019] Lee, N., Ajanthan, T., Torr, P.H.: Snip: Single-shot network pruning based on connection sensitivity. In: International Conference on Learning Representations (ICLR) (2019)
  * Dasgupta and Hsu [2007] Dasgupta, S., Hsu, D.: On-line estimation with the multivariate gaussian distribution. In: International Conference on Computational Learning Theory, pp. 278–292 (2007). Springer
  * Rajasegaran et al. [2020] Rajasegaran, J., Khan, S., Hayat, M., Khan, F.S., Shah, M.: itaml: An incremental task-agnostic meta-learning approach. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 13588–13597 (2020)
  * Salehi et al. [2017] Salehi, S.S.M., Erdogmus, D., Gholipour, A.: Tversky loss function for image segmentation using 3d fully convolutional deep networks. In: International Workshop on Machine Learning in Medical Imaging, pp. 379–387 (2017). Springer
  * Lin et al. [2017] Lin, T.-Y., Goyal, P., Girshick, R., He, K., Dollár, P.: Focal loss for dense object detection. In: Proceedings of the IEEE International Conference on Computer Vision, pp. 2980–2988 (2017)
  * Kingma and Ba [2014] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * Masana et al. [2020] Masana, M., Twardowski, B., Weijer, J.: On class orderings for incremental learning. arXiv preprint arXiv:2007.02145 (2020)


††thanks: This work was supported by the National Natural Science Foundation
of China (Nos. 11805294 and 11975021), the China Postdoctoral Science
Foundation (2018M631013), the Strategic Priority Research Program of Chinese
Academy of Sciences (XDA10010900), the Fundamental Research Funds for the
Central Universities, Sun Yat-sen University (19lgpy268), and in part by the
CAS Center for Excellence in Particle Physics (CCEPP).

# Event vertex and time reconstruction in large-volume liquid scintillator
detectors

Zi-Yuan Li  School of Physics, Sun Yat-Sen University, Guangzhou 510275, China
Yu-Mei Zhang<EMAIL_ADDRESS>Sino-French Institute of Nuclear
Engineering and Technology, Sun Yat-Sen University, Zhuhai 519082, China  Zhen
Qian  School of Physics, Sun Yat-Sen University, Guangzhou 510275, China  Shu
Zhang  School of Physics, Sun Yat-Sen University, Guangzhou 510275, China
Kai-Xuan Huang  School of Physics, Sun Yat-Sen University, Guangzhou 510275,
China  Guo-Fu Cao  Institute of High Energy Physics, Chinese Academy of
Sciences, Beijing 100049, China  Zi-Yan Deng  Institute of High Energy
Physics, Chinese Academy of Sciences, Beijing 100049, China  Gui-Hong Huang
Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049,
China  Wei-Dong Li  Institute of High Energy Physics, Chinese Academy of
Sciences, Beijing 100049, China  Tao Lin  Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, China  Liang-Jian Wen  Institute
of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China
Miao Yu  School of Physics and Technology, Wuhan University, Wuhan 430072,
China  Jia-Heng Zou  Institute of High Energy Physics, Chinese Academy of
Sciences, Beijing 100049, China  Wu-Ming Luo<EMAIL_ADDRESS>Institute of
High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China
Zheng-Yun You<EMAIL_ADDRESS>School of Physics, Sun Yat-Sen
University, Guangzhou 510275, China

###### Abstract

Large-volume liquid scintillator detectors with ultra-low background levels
have been widely used to study neutrino physics and search for dark matter.
Event vertex and event time are not only useful for event selection but also
essential for the reconstruction of event energy. In this study, four event
vertex and event time reconstruction algorithms using charge and time
information collected by photomultiplier tubes were analyzed comprehensively.
The effects of photomultiplier tube properties were also investigated. The
results indicate that the transit time spread is the main effect degrading the
vertex reconstruction, while the effect of dark noise is limited. In addition,
when the event is close to the detector boundary, the charge information
provides better performance for vertex reconstruction than the time
information.

JUNO, Liquid scintillator detector, Neutrino experiment, Vertex
reconstruction, Time reconstruction

##  I Introduction

Liquid scintillators (LSs) have widely been used as detection medium for
neutrinos in experiments such as Kamioka Liquid Scintillator Antineutrino
Detector (KamLAND) kamland , Borexino borexino , Double Chooz doublechooz ,
Daya Bay dyb , and Reactor Experiment for Neutrino Oscillation (RENO) reno .
KamLAND revealed a large mixing angle (LMA) solution for solar neutrino
oscillations. Borexino confirmed the Mikheyev-–Smirnov-–Wolfenstein (MSW) LMA
lmamsw  model in the sub-MeV region for solar neutrino oscillations. Double
Chooz, Daya Bay, and RENO reported nonzero measurements for the mixing angle
$\theta_{13}$. The size of such detectors varies from hundreds to thousands of
cubic meters. Large-volume liquid scintillator detectors are widely used in
the next generation of neutrino experiments, aiming to solve problems such as
neutrinoless double-beta decay (SNO+ at the Sudbury Neutrino Observatory sno+
) and neutrino mass ordering (Jiangmen Underground Neutrino Observatory (JUNO)
junophysics ).

The sensitivity of these experiments is limited by the energy resolution,
detector volume, and detector background. These detectors typically contain a
fiducial volume, where the signal-to-noise ratio is maximal. To distinguish
between events occurring in the fiducial and non-fiducial regions, the event
vertex is reconstructed using the charge and time distribution of photons
collected by the photomultiplier tubes (PMTs). Most importantly, to achieve a
high energy resolution, an accurate vertex is essential to correct for energy
non-uniformity. In addition, the event vertex and event time information can
also be used for particle identification, direction reconstruction, event
classification, and other purposes.

A previous study Liu_2018  investigated vertex reconstruction with time
information in JUNO, without discussing the event time reconstruction, dark
noise effect, and the improvement based on the charge information when the
event is close to the detector boundary. The aim of this study was to analyze
vertex and time reconstruction for point-like events in JUNO under more
realistic conditions. The main contributions of this study are as follows:

  1. $*$

The event time was reconstructed to provide the start time information of the
event, which was important for event alignment, event correlation, etc.

  2. $*$

The dark noise from PMTs was considered and its effect on the vertex
reconstruction was properly controlled.

  3. $*$

Two types of large PMTs were considered and handled separately mainly because
of the difference in transit time spread (TTS).

  4. $*$

An algorithm to provide a more accurate initial vertex value was developed to
improve the performance, especially at the detector boundary region.

  5. $*$

An algorithm employing charge information to reconstruct the event vertex at
the detector boundary was developed.

The remainder of this paper is organized as follows. In Sec. II, a brief
introduction of the JUNO detector and the configurations of the PMTs used in
this study is provided. In Sec. III, the optical processes are described and a
simple optical model is introduced. In Sec. IV, two simple algorithms that can
quickly provide initial values are compared. In Secs. V and VI, two complex
algorithms that provide relatively good performance are introduced in detail.
Finally, a performance summary, the discussion, and the conclusions are
provided in Secs. VII, VIII, and IX, respectively.

##  II The JUNO detector

Figure 1: Schematic of the JUNO detector.

A schematic of the JUNO detector is shown in Fig. 1. The central detector (CD)
is the main part of the JUNO detector assembly, with an acrylic sphere with a
diameter of 35.4 m as inner layer. The acrylic sphere is supported by a
stainless steel latticed shell with a diameter of 40.1 m and filled with
approximately 20,000 tons of LS as target for neutrino detection. The
composition of the LS includes linear alkylbenzene (LAB) as solvent,
2,5-diphenyloxazole (PPO) as fluor, and p-bis-(o-methylstyryl)-benzene (bis-
MSB) as wavelength shifter. To collect photons, the CD is surrounded by
approximately 17,600 20-inch PMTs, including 5000 Hamamatsu R12860 PMTs,
12,600 North Night Vision Technology (NNVT) GDG-6201 PMTs, and approximately
25,600 3-inch XP72B22 PMTs. Around the CD, a 2.5-m-thick water pool is used to
shield external radioactivity from the surrounding rocks and is combined with
approximately 2,000 20-inch GDG-6201 PMTs to serve as a water Cherenkov
detector to veto cosmic ray muons. A top tracker detector, consisting of
plastic scintillators, located above the water pool is used for the
identification and veto of muon tracks. A more detailed description of the
JUNO detector can be found in Refs. junocdr ; junophysics .

The main factors affecting the reconstruction of the event vertex and time
include the TTS and dark noise of the PMTs. In this study, only the 20-inch
PMTs of the CD were used for reconstruction. The number and parameters of the
Hamamatsu and NNVT PMTs are summarized in Table 1 pmtsys . In principle, the
3-inch PMTs with $\sigma\simeq 1.6$ ns TTS could also be included to improve
the reconstruction performance; however, these PMTs 3inchpmt  were not
considered in this study because of their small photon detection coverage
($\sim$3%) gdml ; gdml2 .

Table 1: Number and parameters of the PMTs used in the reconstruction. The TTS and dark noise rate are the mean values of the distribution measured during the mass testing. However, these are not the final values for JUNO. Company | Number | TTS ($\sigma$) | Dark noise rate
---|---|---|---
Hamamatsu | 5,000 | 1.15 ns | 15 kHz
NNVT | 12,600 | 7.65 ns | 32 kHz

##  III Optical processes

When a charged particle deposits energy in the scintillator, the solvent
enters an excited state and transfers energy to the fluor in a non-radiative
manner. Scintillation photons are then emitted along the particle track
through the radiative de-excitation of the excited fluor within a limited
time. The emitted scintillation photons can undergo several different
processes while propagating through a large LS detector. At short wavelengths
($<$ 410 nm), photons are mostly absorbed and then re-emitted at longer
wavelengths, which maximizes the detection efficiency of the PMTs. At long
wavelengths ($>$ 410 nm), photons mainly undergo Rayleigh scattering. A more
detailed study of the wavelength-dependent absorption and re-emission can be
found in Ref. lsopmdl . Additionally, the refractive indices at 420 nm are
1.50 and 1.34 for the LS and water, respectively. The difference in the
refractive indices results in refraction and total reflection at the boundary
of the two media, which affects the time-of-flight of the photons. When using
the time information in the reconstruction, the time-of-flight is crucial, and
it is calculated using the equation

$\mathrm{tof}=\sum_{m}\frac{d_{m}}{v_{m}},$ (1)

where $\mathrm{tof}$, $d_{m}$, and $v_{m}$ are the time-of-flight, optical
path length, and effective light speed, respectively, and $m$ represents
different media, in this case the LS and water, in the JUNO experiment. The
acrylic sphere (thickness of 12 cm) and acrylic cover (thickness of 1 cm) in
front of each PMT were ignored in this study because their refractive indices
were similar to that of the LS and their thickness was small compared to that
of the LS and water.

###  III.1 Optical path length

Figure 2: Event display of the optical path from the event vertex to the PMT
in the JUNO simulation. The red circle ring is the event vertex and the gray
bulbs with blue caps represent the PMTs.

The optical path length can be characterized by the start and end positions in
the detector, which are the vertex of the event and the position of the PMTs,
respectively. The bold cyan curve in Fig. 2 shows a typical example of the
optical path of photons detected by the PMT in the JUNO simulation using the
event display eventdisplay1 ; eventdisplay2 , and further examples are shown
by the thin green curves. There are multiple physically possible paths between
these two positions, each of which has a different optical path length, as
follows:

  1. $*$

owing to absorption and re-emission, the re-emitted photon is not in the same
absorption position and the propagation direction also changes;

  2. $*$

owing to scattering, the photon changes the original direction of the
propagation; and

  3. $*$

owing to refraction and total reflection, the photon does not travel in a
straight line.

As shown in Fig. 2, owing to the various aforementioned optical processes, it
is difficult to predict the actual optical path length for each photon. In
this paper, a simple optical model is proposed, which uses a straight line
connecting the vertex and the PMTs to calculate the optical path length (Fig.
3), and combines with the effective light speed to correct for the time-of-
flight. Using this simple optical model reasonable results can be obtained, as
discussed in Sec. V.

Figure 3: Optical path length from the event vertex to the $i$th PMT. $O$
denotes the center of the detector.

In Fig. 3, {$\vec{r}_{0},t_{0}$} represents the event vertex and start time,
{$\vec{r}_{i},t_{i}$} is the position of the $i$th PMT and the time of the
earliest arriving photon detected by it. The angle between the normal
direction of the $i$th PMT and the vector of the position of the $i$th PMT
pointing to the event vertex is $\theta_{i}$ and
$\alpha_{i}=\arccos(\hat{\vec{r}}_{0}\cdot\hat{\vec{r}}_{i})$. The optical
path length of the photon arriving at the $i$th PMT is
$d_{pathlength,i}=|\vec{r}_{i}-\vec{r}_{0}|=d_{LS,i}+d_{water,i}$ and the
corresponding time-of-flight is $\mathrm{tof}_{i}$. The optical path length in
the LS and water can be calculated by simply solving the trigonometric
equation.

###  III.2 Effective light speed

According to Ref. lsopmdl , the emission spectrum of scintillation photons is
in the range of approximately 300––600 nm. Typically, the group velocity of
the wave packet is used to describe the photon propagation in the medium,
which is given by the equation

$v_{g}(\lambda)=\frac{c}{n(\lambda)-\lambda\frac{\partial{n(\lambda)}}{\partial\lambda}},$
(2)

where $v_{g}$ is the group velocity, $c$ is the speed of light in vacuum, $n$
is the refractive index, and $\lambda$ is the wavelength.

By fitting the Sellmeier equation optics , which describes the dispersion of
the measurement in Refs. rayleighscat  and h2oindex , the refractive index of
the LS and water at different wavelengths is shown in the upper panel of Fig.
4. The group velocity of the LS and water can be calculated using Eq. 2 at
different wavelengths, as shown in the lower panel of Fig. 4.

Figure 4: Dependence of the refractive index (upper panel) and group velocity
(lower panel) on the wavelength in the LS and water.

The propagation speed of photons in water ($v_{water}$) was determined as the
average speed weighted by the probability density function of the photon
wavelength, which was obtained from a Monte Carlo (MC) simulation. As the
absorption and re-emission change the initial wavelength, determining the
propagation speed of photons in the LS ($v_{LS}$) is more complicated. To
consider all wavelength-dependent effects that affect the propagation speed of
photons, the effective light speed $v_{eff}$ is introduced. In addition,
$v_{eff}$ also mitigates the effects by the simplified optical model, which,
for example, ignores the refraction at the interface between the LS and water,
as well as the change in the optical path length due to Rayleigh scattering.
The exact value for $v_{eff}$ can be determined using a data-driven method
based on the calibration data as follows: place $\gamma$ sources along the
Z-axis, use $v_{LS}$ at 420 nm as the initial value of $v_{eff}$ in the
reconstruction algorithm and then, calibrate $v_{eff}$ such that the source
positions can be appropriately reconstructed. As no calibration data was
available for JUNO, in this study, simulated calibration data were used, and
the optimized values for the effective refractive index (c/$v_{eff}$) were
1.546 in the LS and 1.373 in water. In the future, the same method can be
applied to the experimental calibration data.

##  IV Initial value for vertex and time

The TMinuit package tminuit  was used for the minimization procedure in the
time likelihood and in the charge likelihood algorithm introduced in Secs. V
and VI. When there are multiple local minima in the parameter space, an
inaccurate initial value results in local instead of global minima, resulting
in a lower reconstruction efficiency. For detectors such as JUNO, the initial
value needs to be treated carefully because of the total reflection, as
discussed in the following subsections.

###  IV.1 Charge-based algorithm

The charge-based algorithm is essentially based on the charge-weighted average
of the positions of the PMTs in an event, and the event vertex can be
determined using the equation

$\vec{r}_{0}=a\cdot\frac{\sum_{i}{q_{i}\cdot\vec{r}_{i}}}{\sum_{i}{q_{i}}},$
(3)

where $q_{i}$ is the charge of the pulses detected by the $i$th PMT and
$\vec{r}_{0}$ and $\vec{r}_{i}$ are defined in Fig. 3. A scale factor $a$ is
introduced because the charge-based algorithm is inherently biased and an
ideal point-like event in a spherical detector is covered by a uniform
photocathode. Even if all propagation-related effects, such as absorption and
scattering are ignored, the result of a simple integral of the intersections
of all photons with the sphere surface shows that the reconstructed position
of the event is 2/3 of the true position. The value of $a$ can be tuned based
on the calibration data along the $Z$-axis. In this study, $a=1.3$ was used,
which was sufficient to provide an initial estimate for the event vertex.

Figure 5: Heatmap of $R_{rec}$ (upper panel) and $R_{rec}-R_{true}$ (lower
panel) as a function of $R_{true}$ for 4-MeV $e^{+}$ uniformly distributed in
space calculated by the charge-based algorithm.

As can be seen in Fig. 5, even with the scale factor, owing to total
reflection, the reconstructed vertex deviates up to 3 m near the detector
boundary. According to Ref. optics , total reflection occurs only when the
event vertex is located at an $R$ larger than
$R_{LS}\cdot{n_{water}/n_{LS}}\approx 15.9$ m, where $R_{LS}$ is the radius of
the acrylic sphere, $n_{LS}$ and $n_{water}$ are the refractive indices in the
LS and water, respectively. The total reflection region is defined as $R>15.9$
m while $R<15.9$ m is the central region. If the result from the charge-based
algorithm is used as the initial value for the time likelihood algorithm,
approximately 18% of events is reconstructed at a local minimum position. In
addition, it should be noted that the charge-based algorithm is not able to
provide an initial value for the event generation time $t_{0}$. Therefore, a
fast time-based algorithm needs to be introduced, which can provide more
accurate initial values.

###  IV.2 Time-based algorithm

The time-based algorithm uses the distribution of the time-of-flight
correction time $\Delta{t}$ (defined in Eq. 4) of an event to reconstruct its
vertex and $t_{0}$. In practice, the algorithm finds the reconstructed vertex
and $t_{0}$ using the following iterations:

  1. 1.

Apply the charge-based algorithm to obtain the initial vertex.

  2. 2.

Calculate time-of-flight correction time $\Delta{t}$ for the $i$th PMT as

$\Delta{t}_{i}(j)=t_{i}-\mathrm{tof}_{i}(j),$ (4)

where $j$ is the iteration step and $t_{i}$, $\mathrm{tof}_{i}$ are defined in
Fig. 3. Plot the $\Delta{t}$ distribution for all triggered PMTs, and label
the peak position as $\Delta{t}^{peak}$.

  3. 3.

Calculate the correction vector $\vec{\delta}[\vec{r}(j)]$ as

$\vec{\delta}[\vec{r}(j)]=\frac{\sum_{i}(\frac{\Delta{t}_{i}(j)-\Delta{t}^{peak}(j)}{\mathrm{tof}_{i}(j)})\cdot(\vec{r}_{0}(j)-\vec{r}_{i})}{N^{peak}(j)},$
(5)

where $\vec{r}_{0}$, and $\vec{r}_{i}$ are defined in Fig. 3. To minimize the
effect of scattering, reflection, and dark noise on the bias of the
reconstructed vertex, only the pulses appearing in the
$(-10\leavevmode\nobreak\ \rm ns,+5\leavevmode\nobreak\ \rm ns)$ window around
$\Delta{t}^{peak}$ are included. The time cut also suppresses the effect of
the late scintillation photons. The number of triggered PMTs in the window is
$N^{peak}$.

  4. 4.

If $\vec{\delta}[\vec{r}(j)]<1\leavevmode\nobreak\ \rm mm$ or $j=100$, stop
the iteration; otherwise, update the vertex with
$\vec{r}_{0}(j+1)=\vec{r}_{0}(j)+\vec{\delta}[\vec{r}(j)]$ and go to step 2 to
start a new round of iteration.

The distribution of $\Delta{t}$ at different iteration steps is shown in Fig.
6. At the beginning of the iteration, the $\Delta{t}$ distribution is wide
because the initial vertex is far from the true vertex. As the number of
iterations increases, the $\Delta{t}$ distribution becomes more concentrated.
Finally, when the requirement in step 4 is satisfied, the iteration stops. In
the final step, $\vec{r}_{0}$ is the reconstructed vertex and
$\Delta{t}^{peak}$ is the reconstructed time $t_{0}$.

Figure 6: $\Delta{t}$ distribution at different iteration steps $j$.

After the time-of-flight correction, the $\Delta{t}$ distribution is
independent of the event vertex. However, because the earliest arrival time is
used, according to the first-order statistic, as discussed in Ref. fos ; fos2
; fos3 , $t_{i}$ is related to the number of photoelectrons $N^{i}_{pe}$
detected by $i$th PMT. To reduce the bias of the vertex reconstruction, the
following form of the time–$N_{pe}$ correction is applied, and in Eq. 4
$t_{i}$ is replaced by $t^{\prime}_{i}$:

$t^{\prime}_{i}=t_{i}-p0/\sqrt{N^{i}_{pe}}-p1-p2/N^{i}_{pe}.$ (6)

The parameters $(p0,p1,p2)$ with the corresponding values of (9.42, 0.74,
$-$4.60) for Hamamatsu PMTs and (41.31, $-$12.04, $-$20.02) for NNVT PMTs were
found to minimize the bias and energy dependence of the reconstruction in this
study. The difference in the parameters is mainly due to the difference in the
TTSs of the PMTs. Following the correction, the times of different PMTs with
different values of $N_{pe}$ are aligned.

Figure 7: Heatmap of $R_{rec}$ (upper panel) and $R_{rec}-R_{true}$ (lower
panel) as a function of $R_{true}$ for 4-MeV $e^{+}$ uniformly distributed in
space calculated by the time-based algorithm.

As shown in Fig. 7, the time-based algorithm provided a more accurate
reconstructed vertex than the charge-based algorithm (Fig. 5). In addition,
after the time–$N_{pe}$ correction, the reconstruction shows no obvious bias
within the entire detector, even in the total reflection region. The
reconstructed result was used as the initial value for the time likelihood
algorithm.

##  V Time likelihood algorithm

###  V.1 Principle of the algorithm

The time likelihood algorithm uses the scintillator response function to
reconstruct the event vertex. The variable residual time
$t_{res}(\vec{r}_{0},t_{0})$ for the $i$th PMT can be described as

$t^{i}_{res}(\vec{r}_{0},t_{0})=t_{i}-\mathrm{tof}_{i}-t_{0},$ (7)

where $t^{i}_{res}$ is the residual time of the $i$th PMT and $\vec{r}_{0}$,
$t_{0}$, $t_{i}$, and $\mathrm{tof}_{i}$ are defined in Fig. 3.

The scintillator response function mainly consists of the emission time
profile of the scintillation photons and the TTS and the dark noise of PMTs.
In principle, the additional delays introduced by the absorption, re-emission,
scattering, and total reflection of the photon arriving to the PMT depend on
the distance between the emission position and the individual PMTs. However,
the differences are only noticeable for the late arrival hits, which are
largely suppressed by the requirement for the earliest arriving photons in the
time likelihood algorithm. Therefore, in the first-order approximation, the
scintillator response function can be considered to be the same for all
positions inside the scintillator. The scintillator response function can be
described as follows.

As described in Sec. III, when a charged particle interacts with a
scintillator molecule, the molecule is excited, then de-excites, and emits
photons. Typically, the scintillator has more than one component; thus, the
emission time profile of the scintillation photons, $f(t_{res})$, can be
described as

$f(t_{res})=\sum_{k}\frac{\rho_{k}}{\tau_{k}}e^{\frac{-t_{res}}{\tau_{k}}},\sum_{k}\rho_{k}=1,$
(8)

where each ${k}$ component is characterized by its decay time $\tau_{k}$ and
intensity $\rho_{k}$. The different components result from the different
excited states of the scintillator molecules.

To consider the spread in the arrival time of photons at the PMTs, a
convolution with a Gaussian function is applied, given by

$g(t_{res})=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(t_{res}-\nu)^{2}}{2\sigma^{2}}}\cdot{f(t_{res})}.$
(9)

where $\sigma$ is the TTS of PMTs and $\nu$ is the average transit time.

The dark noise, which occurs without incident photons in the PMTs, is not
correlated with any physical event. The fraction of the dark noise in the
total number of photoelectrons $\varepsilon_{dn}$ can be calculated based on
the data acquisition (DAQ) windows, dark noise rate, and light yield of the
LS. The probability of dark noise $\varepsilon(t_{res})$ is constant over
time, where $\int_{DAQ}{\varepsilon(t_{res})dt_{res}}=\varepsilon_{dn}$. By
adding $\varepsilon(t_{res})$ to $g(t_{res})$ and renormalizing its integral
to 1, the probability density function (PDF) of the scintillator response
function can be written as

$p(t_{res})=(1-\varepsilon_{dn})\cdot{g(t_{res})}+\varepsilon(t_{res}).$ (10)

The distribution of the residual time $t_{res}$ of an event for a hypothetical
vertex can be compared with $p(t_{res})$. The best-fitting vertex and $t_{0}$
are chosen by minimizing the negative log-likelihood function

$\mathcal{L}{(\vec{r}_{0},t_{0})}=-\ln(\prod_{i}p(t^{i}_{res})).$ (11)

The parameters in Eq. 10 can be measured experimentally tres1 ; tres2 ; tres3
; tres4 . In this work, the PDF from the MC simulation for the methodology
study was employed.

###  V.2 Probability density function

Figure 8: PDF of the scintillator response function for PMTs detecting
different numbers of photoelectrons. The upper panel shows the response
function for Hamamatsu, the lower panel for the NNVT PMTs.

Figure 9: Bias of the reconstructed $R$ (left panel), $\theta$ (middle panel),
and $\phi$ (right panel) for different energies calculated by the time
likelihood algorithm.

Figure 10: Resolution of the reconstructed $R$ (left panel), $\theta$ (middle
panel), and $\phi$ (right panel) as a function of energy calculated by the
time likelihood algorithm. Different colors represent different PMT
configurations.

The PDF of the scintillator response function for PMTs detecting a single
photoelectron was obtained from the MC simulation, using a 4.4-MeV $\gamma$
source located at the center of the detector, such that the distance to all
PMTs is the same. For PMTs detecting multiple photoelectrons, the time of the
earliest arriving photon is biased toward an earlier time. Therefore, the PDF
need to be modified according to the first-order statistic of $p(t_{res})$ or
the so-called first photoelectron timing technique fos ; fos2 ; fos3  as

$p_{N_{pe}}(t_{res})=N_{pe}p(t_{res})(\int_{t_{res}}^{\infty}p(x)dx)^{N_{pe}-1},$
(12)

where $p_{N_{pe}}(t_{res})$ is the PDF of the scintillator response function
when the PMTs detect $N_{pe}$ hits.

The PDF of two kinds of PMTs is shown in Fig. 8: the upper panel is for
Hamamatsu while the lower panel is for NNVT PMTs. As the PDF is affected by
the time resolution of the PMTs, the PDF of the NNVT is wider because of its
inferior TTS. The inset in the lower panel shows the PDF on a logarithmic
scale, and the time constant contribution of the dark noise
$\varepsilon(t_{res})$ is clearly visible.

###  V.3 Reconstruction performance

The reconstructed vertex was compared with the true vertex in spherical
coordinates ($R,\theta,\phi$) for the MC $e^{+}$ samples and fitted with a
Gaussian function to analyze the bias and resolution. The bias of the
reconstruction is shown in Fig. 9, where different colors represent events
with different energies. As can be seen in the left panel of Fig. 9, the
reconstructed $R$ is consistent with the true value in the central region,
while an energy-dependent bias behavior is noticeable near the detector
boundary. Given its regular bias behavior, the bias can be corrected with an
energy-dependent correction. Moreover, although the reconstructed $R$ is
biased, there is no bias in $\theta$ and $\phi$, as shown in the middle and
right panels of Fig. 9, respectively.

The spatial resolution of the vertex reconstruction as a function of energy is
shown in Fig. 10. The $R$ bias was corrected before the analysis of the
resolution. To study the individual effect of the TTS and dark noise on the
vertex reconstruction, different MC samples were produced with and without
these effects. The vertex reconstruction results are shown in Fig. 10. The
magenta circles represent the default PMT configuration, as described in Sec.
II. The red triangles represent an ideal configuration, which assumes perfect
PMTs without the effects of the TTS and dark noise. The black squares
represent the configuration of PMTs including only the dark noise effect,
while the blue inverted triangles represent the PMT configuration including
only the TTS effect. The exact values of the vertex resolution at 1.022 MeV
and 10.022 MeV are summarized in Tables 2 and 3, respectively. The energy
$E_{true}$ includes the energy of the annihilation gamma rays. The light yield
was approximately 1300 detected $N_{pe}$ per 1 MeV of deposited energy in
JUNO, and the energy nonlinearities on the light yield were ignored in the
approximation. As can be seen in Tables 2 and 3, the dark noise has no effect
at high energy and its effect at low energy is also highly limited. The
largest effect results from the TTS in the time likelihood algorithm. The
energy-dependent vertex resolution is approximately proportional to
$1/\sqrt{N_{pe}}$ fos2 .

Table 2: Vertex resolution for different PMT configurations at 1.022 MeV (detection of $\sim$1328 $N_{pe}$ in total, corresponding to $\sim$370 $N_{pe}$ detected by Hamamatsu PMTs). PMT configuration |  $R$ (mm) |  $\theta$ (degrees) |  $\phi$ (degrees)
---|---|---|---
Ideal | 60 | 0.25 | 0.31
With dark noise only | 62 | 0.27 | 0.34
With TTS only | 89 | 0.37 | 0.44
With TTS and dark noise | 103 | 0.40 | 0.47
With TTS and dark noise (Hamamatsu PMTs only) | 105 | 0.42 | 0.49
Table 3: Vertex resolution for different PMT configurations at 10.022 MeV (detection of $\sim$13280 $N_{pe}$ in total, corresponding to $\sim$ 3700 $N_{pe}$ detected by Hamamatsu PMTs). PMT configuration |  $R$ (mm) |  $\theta$ (degrees) |  $\phi$ (degrees)
---|---|---|---
Ideal | 19 | 0.08 | 0.11
With dark noise only | 19 | 0.08 | 0.11
With TTS only | 31 | 0.13 | 0.16
With TTS and dark noise | 31 | 0.13 | 0.16
With TTS and dark noise (Hamamatsu PMTs only) | 32 | 0.14 | 0.17

Owing to the low time resolution of the NNVT PMTs, in Fig. 10 only the
reconstruction using Hamamatsu PMTs is shown (green circles). In this study,
we found that the vertex resolution with Hamamatsu PMTs was similar to that of
using all PMTs. The reconstruction speed was 3.5 times faster, because the
fraction of the Hamamatsu PMTs was approximately 28% of all PMTs in the CD.

Figure 11: Reconstructed event time $t_{0}$ at different energies.

The reconstructed event time $t_{0}$ is shown in Fig. 11. The effect of
$t_{0}$ is essentially a global shift of an event to match the scintillator
response function PDF; in reality, $t_{0}$ is also affected by the trigger
time and the time delay from the cable. The absolute value of $t_{0}$ can be
neglected; only the relative difference of different events is important for
the alignment of events. The small bump near $-$1.6 ns is correlated with the
$R$ bias, and the long tail on the right side results from positronium
formation. The variation in the reconstructed $t_{0}$ is within a few
nanoseconds.

##  VI Total reflection region calculated by the charge likelihood algorithm

The time likelihood method described in Sec. V introduces a bias in the $R$
direction when the reconstructing events are close to the acrylic sphere. As
mentioned in Ref. Smirnov_2003 , using a charge signal with the maximum
likelihood method can provide better spatial resolution than the time
likelihood algorithm when an event occurs near the detector boundary. In this
section, we discuss the charge likelihood algorithm to reconstruct the event
vertex in the total reflection region only, while the reconstruction result in
the central region is omitted.

The charge likelihood algorithm is based on the distribution of the number of
photoelectrons in each PMT. With the mean expected number of photoelectrons
$\mu(\vec{r_{0}},E)$ detected by each PMT at a given vertex and energy, the
probability of observing $N_{pe}$ on a PMT follows a Poisson distribution.
Furthermore,

  1. $*$

Probability for the $j$th PMT with no hits:
$P_{nohit}^{j}(\vec{r_{0}},E)=e^{-\mu_{j}}$,

  2. $*$

Probability for the $i$th PMT with $N^{i}_{pe}$ hits:
$P_{hit}^{i}(\vec{r_{0}},E)=\frac{\mu_{i}^{N^{i}_{pe}}e^{-\mu_{i}}}{N^{i}_{pe}!}$.

Therefore, the probability of observing a hit pattern for an event can be
written as

$p(\vec{r_{0}},E)=\prod_{j}{P_{nohit}^{j}(\vec{r_{0}},E)}\cdot{\prod_{i}{P_{hit}^{i}(\vec{r_{0}},E)}}.$
(13)

The best-fit values of $\vec{r_{0}}$ and $E$ can be obtained by minimizing the
negative log-likelihood function

$\mathcal{L}{(\vec{r_{0}},E)}=-\ln(p(\vec{r_{0}},E)).$ (14)

In principle, $\mu(\vec{r_{0}},E)$ can be expressed by the equation

$\mu_{i}(\vec{r_{0}},E)=Y\cdot\frac{\Omega(\vec{r_{0}},r_{i})}{4\pi}\cdot\varepsilon_{i}\cdot{f(\theta_{i})}\cdot{e^{-\sum_{m}{\frac{d_{m}}{\zeta_{m}}}}}\cdot{E}+\delta_{i},$
(15)

where $Y$ is the energy scale factor, $\Omega(\vec{r_{0}},r_{i})$ is the solid
angle of the $i$th PMT, $\varepsilon_{i}$ is the detection efficiency of the
$i$th PMT, $f(\theta_{i})$ is the angular response of the $i$th PMT,
$\theta_{i}$ is defined in Fig. 3, $\zeta_{m}$ is the attenuation length
attleng  in materials, and $\delta_{i}$ is the expected number of dark noise.
This equation is based on the assumption that the scintillation light yield is
linearly proportional to the energy.

Figure 12: Mean expected number of photoelectron distribution as a function of
radius $R$ and angle $\alpha$. This map is obtained by placing gamma sources
at 29 specific positions along the Z-axis, which can be performed using a
calibration procedure Wu_2019 .

However, Eq. 15 cannot describe properly the contribution of the indirect
light, the effect of light shadows because of the geometric structure, and the
effect of the total reflection. Another solution is to use the model-
independent method described in Ref. Wu_2019 : the mean expected number of
photoelectrons can be obtained by placing gamma sources at 29 specific
positions along the Z-axis, which can be performed using a calibration
procedure calibsys . In this study, different from Ref. Wu_2019 , we focused
on the performance of the vertex reconstruction. The mean expected number of
the photoelectron distributions as a function of radius $R$ and angle $\alpha$
is shown in Fig. 12, and the definition of angle $\alpha$ is shown in Fig. 3.

The mean expected number of photoelectrons $\mu$ obtained from Fig. 12 was
used to calculate the hit probability. Instead of reconstructing
($R,\theta,\phi$) at the same time, $\theta$ and $\phi$ were fixed at the
reconstructed values provided by the time likelihood algorithm, and only the
event radius $R$ was reconstructed using the charge likelihood algorithm.
Therefore, the probability in Eq. 13 can be rewritten as

$p(R,E)=\prod_{j}{P_{nohit}^{j}(R,E)}\cdot{\prod_{i}{P_{hit}^{i}(R,E)}}.$ (16)

The reconstruction performance, focusing on the total reflection region, is
shown in Figs. 13 and 14. In the total reflection region, the mean value of
the reconstructed $R$ was consistent with the true $R$, and the resolutions in
the $R$ direction were 81 mm at 1.022 MeV and 30 mm at 10.022 MeV.

Figure 13: Bias of the reconstructed $R$ in the total reflection region at
different energies calculated by the charge likelihood algorithm.

Figure 14: Resolution of reconstructed $R$ as a function of energy calculated
by the time likelihood and charge likelihood algorithms
($R^{3}\leavevmode\nobreak\ >\leavevmode\nobreak\ 4000\leavevmode\nobreak\
\rm{m}^{3}$)

.

As the charge distribution provides good radial discrimination ability, this
algorithm can provide better resolution and a significantly smaller bias
compared with those of the time likelihood algorithm in the total reflection
region.

##  VII Performance summary

The execution time of the reconstruction for each event was tested on a
computing cluster with Intel Xeon Gold 6238R CPUs (2.2 GHz), as shown in Fig.
15. The execution time of the charge-based algorithm was in the order of
$O(10^{-4})$ s per event, which cannot be presented in the figure. The
execution time of the time-based and the time likelihood algorithm was
proportional to the event energy and could be reduced by using only the
Hamamatsu PMTs for the reconstruction. The execution time of the charge
likelihood algorithm was independent of the event energy.

Figure 15: Execution time for the reconstruction for different algorithms.

The resolutions of the four algorithms in the $R$ direction are shown in Fig.
16. Owing to the large bias of the charge-based algorithm, a correction to
remove the position-dependent bias was applied before the analysis of the
resolution.

Figure 16: Resolution of the reconstructed $R$ as a function of energy for
different algorithms.

The charged-based algorithm is suitable for online reconstruction tasks that
require high speed but do not require high resolution. The time-based
algorithm does not rely on MC; it can be used as a data-driven reconstruction
method. The time likelihood and charge likelihood algorithms are relatively
accurate and each has its own advantages in a specific detector region.

##  VIII Discussion

The vertex resolutions of KamLAND and Borexino are approximately 12 cm and 10
cm at 1 MeV, respectively, and for JUNO it is approximately 10.5 cm. The
diameter of JUNO (35.4 m) is several times larger than that of KamLAND (13 m)
and Borexino (8.5 m). Despite its larger size, JUNO is still able to achieve a
similar vertex resolution based only on the PMT time information. In this
study, various effects on the vertex reconstruction for JUNO were
comprehensively analyzed. As expected, the TTS of the PMT is the dominant
factor. The vertex reconstruction capability of JUNO is mainly based on the
Hamamatsu PMTs. Although the number of NNVT PMTs is more than twice, their
time information is not useful in the vertex reconstruction because of their
significantly inferior TTS. After considering the light yield of the LS and
the PMT coverage, the number of photons detected by the Hamamatsu PMTs is in
the same range as those of the three detectors mentioned above. This provides
an explanation of their similar vertex resolutions based only on the PMT time
information. To fully exploit the large PMT coverage and large number of PMTs
in JUNO, the charge information of PMTs also need to be utilized in addition
to the time information to constrain the event vertex, especially near the
detector boundary. At the same time, the effect of the dark noise of the PMT
can be mitigated with appropriate treatment. A more accurate initial value
also improves the performance of the vertex reconstruction. In addition to the
event vertex, the event time can also be reconstructed simultaneously, which
is a useful variable for downstream analyses.

In LS detectors, in addition to the scintillation photons there are also
Cherenkov photons, whose effects need to be studied in the future. The $R$
bias near the detector edge in the current results also indicates that a more
accurate PDF of the scintillator response function is needed to include its
dependence on the position as well as the particle type. All vertex
reconstruction methods based on PMT time information use the time of the first
arrival photons only; in principle, later photons might be useful as well.
Therefore, novel methods also need to be explored. Our preliminary studies on
algorithms based on machine learning qianz_2021  showed comparable vertex
reconstruction performance, and they need to be further investigated

##  IX Conclusion

In this study, four algorithms for the reconstruction of the event vertex and
event time were investigated in detail and verified using MC samples generated
by the offline software of JUNO. Considering the TTS and dark noise effects
from the PMTs, a vertex resolution of 10 cm or higher can be achieved in the
energy range of reactor neutrinos. The TTS has a major effect on the vertex
resolution, whereas the effect of dark noise is limited. Near the detector
boundary, charge information can constrain the event vertex better compared
with time information. The algorithms discussed in this paper are applicable
to current and future experiments using similar detection techniques.

## References

  * (1) KamLAND Collaboration, First Results from KamLAND: Evidence for Reactor Antineutrino Disappearance. Phys. Rev. Lett. 90 (2003) 021802. https://doi.org/10.1103/PhysRevLett.90.021802
  * (2) BOREXINO Collaboration, Neutrinos from the primary proton–proton fusion process in the sun. Nature 512 (2014) 383–386. https://doi.org/10.1038/nature13702
  * (3) DOUBLE CHOOZ Collaboration, Indication of reactor ${\overline{\nu}}_{e}$ disappearance in the Double CHOOZ experiment. Phys. Rev. Lett. 108 (2012) 131801. https://doi.org/10.1103/PhysRevLett.108.131801
  * (4) DAYA BAY Collaboration, Observation of electron-antineutrino disappearance at Daya Bay. Phys. Rev. Lett. 108 (2012) 171803. https://doi.org/10.1103/PhysRevLett.108.171803
  * (5) RENO Collaboration, Observation of reactor electron antineutrinos disappearance in the RENO experiment. Phys. Rev. Lett. 108 (2012) 191802. https://doi.org/10.1103/PhysRevLett.108.191802
  * (6) S.T. Petcov, M. Piai, The LMA MSW solution of the solar neutrino problem, inverted neutrino mass hierarchy and reactor neutrino experiments. Phys. Lett. B 533 (2002) 94–106. https://doi.org/10.1016/S0370-2693(02)01591-5
  * (7) SNO+ Collaboration, Current status and future prospects of the SNO+ experiment. Adv. High Energy Phys. 2016 (2016) 6194250. https://doi.org/10.1155/2016/6194250
  * (8) JUNO Collaboration, Neutrino physics with JUNO. J. Phys. G 43 (2016) 030401. https://doi.org/10.1088/0954-3899/43/3/030401
  * (9) Q. Liu, M. He, X. Ding et al., A vertex reconstruction algorithm in the central detector of JUNO. JINST 13 (09) (2018) T09005. https://doi.org/10.1088/1748-0221/13/09/t09005
  * (10) T. Lin, J.H. Zou, W.D. Li et al., The application of SNiPER to the JUNO simulation. J. Phys. Conf. Ser. 898 (2017) 042029. https://doi.org/10.1088/1742-6596/898/4/042029
  * (11) JUNO Collaboration, JUNO Conceptual Design Report. arXiv:1508.07166 [physics.ins-det]
  * (12) Z.M. Wang, JUNO PMT system and prototyping. J. Phys. Conf. Ser. 888 (2017) 012052. https://doi.org/10.1088/1742-6596/888/1/012052
  * (13) G. Wang, Y.K. Heng, X.H. Li et al., Effect of divider current on 7.62cm photomultiplier tube performance. Nucl. Tech. 2018, 41(8): 80402-080402. https://doi.org/10.11889/j.0253-3219.2018.hjs.41.080402
  * (14) K.J. Li, Z.Y. You, Y.M. Zhang et al., GDML based geometry management system for offline software in JUNO. Nucl. Instrum. Meth. A 908 (2018) 43–48. https://doi.org/10.1016/j.nima.2018.08.008
  * (15) S. Zhang, J.S. Li, Y.J. Su et al., A method for sharing dynamic geometry information in studies on liquid-based detectors. Nucl. Sci. Tech. 32, 21 (2021). https://doi.org/10.1007/s41365-021-00852-8
  * (16) Y. Zhang, Z.Y. Yu, X.Y. Li et al., A complete optical model for liquid-scintillator detectors. Nucl. Instrum. Meth. A 967 (2020) 163860. https://doi.org/10.1016/j.nima.2020.163860
  * (17) Z.Y. You, K. Li, Y. Zhang et al., A ROOT based event display software for JUNO, JINST 13 (2018) T02002. https://doi.org/10.1088/1748-0221/13/02/t02002
  * (18) J. Zhu, Z. You, Y. Zhang et al., A method of detector and event visualization with unity in JUNO. JINST 14 (01) (2019) T01007. https://doi.org/10.1088/1748-0221/14/01/t01007
  * (19) M. Born et al., Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (7th ed.). Cambridge: Cambridge University Press, (1999).
  * (20) X. Zhou, Q. Liu, M. Wurm et al., Rayleigh scattering of linear alkylbenzene in large liquid scintillator detectors. Rev. Sci. Instrum. 86 (04 2015). https://doi.org/10.1063/1.4927458
  * (21) A.N. Bashkatov, E.A. Genina, Water refractive index in dependence on temperature and wavelength: a simple approximation. International Society for Optics and Photonics, SPIE, 2003, pp. 393 – 395. https://doi.org/10.1117/12.518857
  * (22) F. James, MINUIT Function Minimization and Error Analysis: Reference Manual Version 94.1 (1994).
  * (23) G. Ranucci, An analytical approach to the evaluation of the pulse shape discrimination properties of scintillators. Nucl. Instrum. Meth. A 354 (2) (1995) 389 – 399. https://doi.org/10.1016/0168-9002(94)00886-8
  * (24) C. Galbiati, K. McCarty, Time and space reconstruction in optical, non-imaging, scintillator-based particle detectors. Nucl. Instrum. Meth. A 568 (2) (2006) 700 – 709. https://doi.org/10.1016/j.nima.2006.07.058
  * (25) M. Moszynski, B. Bengtson, Status of timing with plastic scintillation detectors. Nucl. Instrum. Meth. A 158 (1979) 1 – 31. https://doi.org/10.1016/S0029-554X(79)90170-8
  * (26) W.L. Zhong, Z.H. Li, C.G. Yang et al., Measurement of decay time of liquid scintillator. Nucl. Instrum. Meth. A 587 (2) (2008) 300 – 303. https://doi.org/10.1016/j.nima.2008.01.077
  * (27) T.M. Undagoitia, F. von Feilitzsch, L. Oberauer et al., Fluorescence decay-time constants in organic liquid scintillators. Rev. Sci. Instrum. 80 (4) (2009) 043301. https://doi.org/10.1063/1.3112609
  * (28) H.M. O’Keeffe, E. O’Sullivan, M.C. Chen, Scintillation decay time and pulse shape discrimination in oxygenated and deoxygenated solutions of linear alkylbenzene for the SNO+ experiment. Nucl. Instrum. Meth. A 640 (1) (2011) 119 – 122. https://doi.org/10.1016/j.nima.2011.03.027
  * (29) H. Takiya, K. Abe, K. Hiraide et al., A measurement of the time profile of scintillation induced by low energy gamma-rays in liquid xenon with theXMASS-I detector. Nucl. Instrum. Meth. A 834 (2016) 192 – 96. https://doi.org/10.1016/j.nima.2016.08.014
  * (30) O.J. Smirnov, Energy and Spatial Resolution of a Large-Volume Liquid-Scintillator Detector. Instrum. Exp. Tech. 46 (3) (2003) 327 – 344. https://doi.org/10.1023/a:1024458203966
  * (31) R. Zhang, D.W. Cao, C.W. Loh et al. Using monochromatic light to measure attenuation length of liquid scintillator solvent LAB. Nucl. Sci. Tech. 30, 30 (2019). https://doi.org/10.1007/s41365-019-0542-1
  * (32) W.J. Wu, M. He, X. Zhou et al., A new method of energy reconstruction for large spherical liquid scintillator detectors. JINST 14 (03) (2019) P03009. https://doi.org/10.1088/1748-0221/14/03/p03009
  * (33) G.L. Zhu, J.L. Liu, Q. Wang et al., Ultrasonic positioning system for the calibration of central detector. Nucl. Sci. Tech. 30, 5 (2019). https://doi.org/10.1007/s41365-018-0530-x
  * (34) Z. Qian, V. Belavin, V. Bokov et al., Vertex and Energy Reconstruction in JUNO with Machine Learning Methods. arXiv:2101.04839 [physics.ins-det]