SECLEDS: Sequence Clustering in Evolving Data Streams via Multiple Medoids and
Medoid Voting 11institutetext: Delft University of Technology, The Netherlands
11email<EMAIL_ADDRESS>
###### Abstract

Sequence clustering in a streaming environment is challenging because it is
computationally expensive, and the sequences may evolve over time. K-medoids
or Partitioning Around Medoids (PAM) is commonly used to cluster sequences
since it supports alignment-based distances, and the k-centers being actual
data items helps with cluster interpretability. However, offline k-medoids has
no support for concept drift, while also being prohibitively expensive for
clustering data streams. We therefore propose SECLEDS, a streaming variant of
the k-medoids algorithm with constant memory footprint. SECLEDS has two unique
properties: i) it uses multiple medoids per cluster, producing stable high-
quality clusters, and ii) it handles concept drift using an intuitive Medoid
Voting scheme for approximating cluster distances. Unlike existing adaptive
algorithms that create new clusters for new concepts, SECLEDS follows a
fundamentally different approach, where the clusters themselves evolve with an
evolving stream. Using real and synthetic datasets, we empirically demonstrate
that SECLEDS produces high-quality clusters regardless of drift, stream size,
data dimensionality, and number of clusters. We compare against three popular
stream and batch clustering algorithms. The state-of-the-art BanditPAM is used
as an offline benchmark. SECLEDS achieves comparable F1 score to BanditPAM
while reducing the number of required distance computations by 83.7%.
Importantly, SECLEDS outperforms all baselines by 138.7% when the stream
contains drift. We also cluster real network traffic, and provide evidence
that SECLEDS can support network bandwidths of up to 1.08 Gbps while using the
(expensive) dynamic time warping distance.

###### Keywords:

Sequence clustering K-medoids Data streams Concept drift Network traffic.

##  1 Introduction

Stream clustering is the problem of clustering a potentially unbounded stream
of items in a single pass, where the items arrive sequentially without any
particular order, _e.g.,_ network traffic, financial transactions, and sensor
data. Stream clustering algorithms must have low memory overhead, be
computationally efficient, and robust to concept drift, _i.e.,_ evolving data
distributions [25]. Maintaining high cluster quality in a fully online setting
is extremely difficult. Therefore, hybrid online-offline algorithms are
popular among existing approaches, _e.g.,_ CluStream [2], StreamKM++ [1],
DenStream [6], BIRCH [31]. These algorithms have an online component that
summarizes the data stream, and an offline component that periodically uses
that information to create the final clusters. There also exist algorithms
that store part of the stream for handling outliers, _e.g.,_ BOCEDS [15], MDSC
[9]. Existing stream clustering algorithms handle concept drift by having
variable number of clusters: they add new clusters for newly observed behavior
and discard clusters that contain too many old data items. This leads to
higher memory requirements for managing buffers and intermediate solutions.
Batch clustering algorithms can also be used in a streaming setting by
considering a batch size of one, _e.g.,_ Minibatch k-means [24]. However, they
start to under-perform when the stream contains drift.

In recent years, sequential data has increasingly become popular because of
the powerful insights that it provides regarding behavior analytics [5],
_e.g.,_ for attacker strategy profiling [21], fraud detection [13], human
activity recognition [7]. Clustering sequences in an offline setting is
challenging in itself because sequences are often out-of-sync, requiring
expensive alignment-based distance measures, which are often not supported by
many clustering algorithms. K-medoids or Partitioning Around Medoids (PAM) has
often been used to cluster sequences because the k-centers are represented by
actual data items, called medoids or prototypes [28, 27]. This has multiple
benefits: i) it makes the cluster interpretation simpler; ii) it enables the
use of non-metric distances such as dynamic time warping (DTW); and iii) it
allows to estimate exact storage requirements based on the k-fixed clusters.
Although the state-of-the-art offline k-medoids algorithms, _i.e.,_ FastPAM1
[23] and BanditPAM [26] have reduced the runtime complexity to
$\mathcal{O}(nlogn)$, they are still not efficient enough to be used in
streaming settings, and the cluster quality will degrade over time as the
stream evolves. To the best of our knowledge, there exists no streaming
version of the k-medoids algorithm that can efficiently cluster sequential
data.

Figure 1: An illustration of SECLEDS’ clusters following an evolving data
stream. The medoids close to recent data gain more votes, while the medoids
with the least votes are replaced with new data items from the stream. In
effect, the k-clusters handle concept drift by capturing different concepts in
the stream at different time steps.

Contributions. In this paper, we propose SECLEDS, a lightweight streaming
version of the k-medoids algorithm with constant memory footprint. SECLEDS has
two unique properties: Firstly, it uses p-medoids per cluster to maintain
stable high-quality clusters. Note the difference from IMMFC [30], which uses
the information of multiple medoids in independent sub-solutions to select the
final medoids. We initialize the p-medoids using a non-uniform sampling
strategy similar to k-means++. Secondly, a Medoid Voting scheme is used to
estimate a cluster’s center of mass. The offline k-medoids has a Swap step
that tests each point in a cluster to determine the next medoid. SECLEDS
cannot do this because it does not store any part of the stream. Instead, it
maintains votes for each medoid that estimate how representative (valuable) it
is given the data seen so far. A user-supplied decay factor enables SECLEDS to
slowly forget the votes regarding past data. The least representative medoids
are then replaced with new data items. This way, rather than creating new
clusters for new concepts, the k-clusters themselves evolve with the data
stream. Figure 1 shows how the clusters follow a data stream as it evolves.
Thus, the k-clusters represent different concepts in the stream at different
time steps. SECLEDS addresses the following real-world constraints:

  1. I.

A runtime efficient medoid-based clustering algorithm with a fixed memory
footprint that can handle high-bandwidth data streams,

  2. II.

An algorithm that produces high-quality clusters in a streaming environment
while being able to deal with concept drift,

  3. III.

Accurate sequence clusters using alignment-based distances, and

  4. IV.

Minimal parameter settings to support ease-of-use.

Empirical results.  We experiment on several real and synthetic data streams
that contain 2D points and univariate sequences. We empirically demonstrate
that SECLEDS produces high-quality clusters regardless of drift, stream size,
and number of clusters. We use the following state-of-the-art and popular
clustering algorithms as baselines: a) Streaming: CluStream, StreamKM++; b)
Batch: Minibatch k-means; c) Offline: BanditPAM. Particularly, BanditPAM is
used as a benchmark for the best achievable clustering on a static dataset.
The results show that i) SECLEDS achieves comparable F1 score to BanditPAM,
while reducing the required number of distance computations by 83.7%; ii)
SECLEDS outperforms all baselines by 138.7% when the stream contains drift;
iii) SECLEDS is faster than BanditPAM and CluStream on most clustering tasks.

We also discuss a use-case from the network security domain where network
traffic is often randomly sampled to keep the storage requirements within a
predefined budget. Consequently, temporal patterns in the network traffic are
lost that could have been useful for downstream tasks, _e.g.,_ behavior
analytics. We propose a smarter sampling technique that uses medoid-based
stream clustering (SECLEDS) to summarize the network traffic: SECLEDS clusters
sequences of network traffic and periodically stores the medoids of each
cluster, thus reducing the storage needs while preserving temporal patterns in
the data. By clustering real-world network traffic, we provide evidence that
SECLEDS (and SECLEDS-dtw) can support network bandwidths of up to 2.79 Gbps
(and 1.08 Gbps), respectively. We release SECLEDS as open-source111SECLEDS:
https://github.com/tudelft-cda-lab/SECLEDS.

##  2 Preliminaries

Stream. Given a sensor that receives an unbounded stream of multi-dimensional
data points $X=\\{x_{1},x_{2},\dots\\}$ with dimensionality $d$, arriving at
time steps $T=\\{t_{1},t_{2},\dots\\}$, a sequential data stream is defined as
$S=\\{s_{1},\dots s_{n},\dots\\}$, where $s_{i}$ is a time window $w$ over $X$
such that $s_{i}=\\{x_{i},x_{i+1},\dots x_{i+w}\\}$, and $y_{i}$ is its
associated class label. Traditional point clustering considers $w=1$, while
for sequence clustering, we consider $w>1$. We use two configurations, _i.e.,_
$d=2,w=1$ (2D point clustering) and $d=1,w=100$ (univariate sequence
clustering). A case of bivariate variable length sequences is given in
appendix.

Concept drift. Real-world data streams often change unexpectedly over time.
This shift alters the statistical properties of their underlying distribution.
In machine learning, this is called concept drift [18, 32]. Concept drift is
typically categorized into four types [17]: (i) Sudden drift where a new
concept arises abruptly; (ii) Gradual drift where an old concept is slowly
replaced by a new one; (iii) Incremental drift where a concept incrementally
turns into another one; and (iv) Recurring concepts are old concepts that
reappear from time to time. Years of research has gone into developing concept
drift detectors that either monitor the underlying data distribution, error
rate or perform hypothesis testing to trigger model retraining [17, 4].
Typical stream clustering algorithms handle concept drift by introducing new
clusters for new concepts, and discarding old irrelevant clusters [14].
Although intuitively appealing, this requires user-supplied parameters that
define what ‘new’ means.

##  3 SECLEDS: Sequence Clustering in Evolving Streams

SECLEDS is a lightweight streaming variant of the classical k-medoids (PAM)
algorithm. To support high bandwidth data streams, SECLEDS does not store any
part of the stream in memory — it receives an item, assigns it to one of the
k-clusters, and then discards it. This way, SECLEDS has a guaranteed constant
memory footprint [I]. However, this requirement cannot be achieved using the
offline Build and Swap steps of PAM. Instead, SECLEDS performs a non-uniform
sampling (similar to k-means++) on an initial batch of the stream to
initialize the medoids. It also makes use of multiple medoids per cluster to
provide a stable cluster definition in a streaming setting where noise and
concept drift are common properties [II].

The efficiency of a chosen distance measure is usually the primary performance
bottleneck in sequence clustering. Thus, minimizing the number of distance
computations is key to scaling SECLEDS to large data streams. This is achieved
by introducing a Medoid Voting scheme whose purpose is twofold: (i) it
determines the center of mass of a cluster, thus is able to estimate how
representative a medoid is given the data seen so far; and (ii) by using this
information, old irrelevant medoids are replaced by new ones that are located
near recent data. Hence, better medoids can be found without having to perform
additional distance computations [I]. This also allows SECLEDS to support
robust but computationally expensive distance measures specifically meant for
sequential data, _e.g.,_ dynamic time warping [III]. Finally, SECLEDS handles
concept drift by regularly forgetting past data and occupying newer
regions/concepts in the data stream. This is achieved by applying exponential
decay $\lambda$ to the medoid votes at each time step [II].

SECLEDS has a modular implementation in Python. The k-clusters, p-medoids per
cluster, and decay rate $\lambda$ are the only three user-supplied parameters
needed for the algorithm, making it useful for exploratory data analysis [IV].
We believe these parameters are easier to tune compared to many radius- or
density-based hyperparameters in existing clustering algorithms, which require
a deeper understanding of the data distribution in advance.

Figure 2: The effect of cluster initialization and concept drift on medoid
selection of a cluster, where $p=5$. (Top): Given a uniformly distributed
stream, every medoid becomes popular at some point. (Middle): For an
incrementally drifted stream, the medoid close to drifted data becomes
popular. (Bottom): For a class-ordered stream, all clusters start from one
class, until one medoid migrates to the correct class. In this case, the
correct class is observed from $t_{1300}$. Medoid-3 migrates first and becomes
popular, while medoid-1 migrates last.

###  3.1 Stable cluster definition via multiple medoids

A new data item $s$ is assigned to a cluster $cid$ with the least average
distance to its medoids. With multiple medoids per cluster, this provides a
robust cluster assignment. Additionally, the medoid voting scheme encourages
the medoids to represent different sections of a class so they can gain votes:
if they are too close together, some of them do not receive votes and get
replaced eventually. It also ensures that outliers are quickly replaced
because of fewer votes.

Concept drift and cluster initialization determine how the medoids behave.
Figure 2 illustrates three scenarios with varying medoid behavior for a single
cluster as a function of votes gained over time: (a) Assuming no concept
drift, when the stream is roughly evenly shuffled, all the medoids receive
uniform votes on average. This is because all medoids are close to parts of
the stream at different time steps. At a specific time step, the medoid close
to the most amount of recent data becomes popular. The top figure shows that
each medoid becomes popular at some point in the stream, indicating that the
medoids represent different sections of the underlying class. (b) When the
stream incrementally drifts, the cluster follows the evolving stream by
replacing the least popular medoids with recent data items from the stream.
Since the new medoids are now closer to new data, they gain more votes and
become popular. This has roughly the same effect as the first case. (c) When
the data arrives one class at a time, all clusters are initialized in a single
class. As data from a new class appears, one medoid from the closest cluster
migrates to it and starts gaining votes. Over time, the older popular medoids
lose their votes because of exponential decay, and eventually migrate to the
new class. This is shown from $t_{1300}$ onward in the bottom figure,
highlighting the importance of multiple medoids in noisy streams.

###  3.2 Center of mass estimation

The voting scheme provides an estimate of a cluster’s center of mass by
assigning more votes to recently observed data in the stream $S$, while
exponential decay helps to forget votes regarding older data. Without decay,
older clusters with popular medoids never evolve. Thus, these properties help
to replace irrelevant medoids, _e.g.,_ those that are located i) close to the
least amount of recent data, or ii) in a region where new data no longer
arrives. Note that we only apply exponential decay to the most recently
updated cluster, so that we do not forget valuable information about other
clusters while the data from this class arrives.

##  4 The SECLEDS Algorithm

SECLEDS has three modules: an initialization module (Init), an assignment
module (Assign), and an update module (Update). The task is to assign each
item in $S$ to one of the k-clusters. SECLEDS maintains and updates a model of
the stream seen so far in the form of k-clusters, $\mathcal{C}=\\{C_{1},\dots
C_{k}\\}$. For clarity, we use $t$ to denote the clusters at time $t$. These
superscripts are removed from Algorithm 1. Each cluster is represented by a
set of p-medoids and their votes, _i.e.,_ for $1\leq i\leq k$,
$C_{i}^{t}=\\{(m^{t}_{i,1},v^{t}_{i,1})\ldots(m^{t}_{i,p},v_{i,p}^{t})\\}$,
where for $1\leq j\leq p$: $m^{t}_{i,j}\in S$ is the $j^{th}$ medoid of the
$i^{th}$ cluster at time $t$ having $v_{i,j}^{t}\in\mathbb{R}$ votes.

Init. A batch $\mathbb{B}$ from the start of $S$ is used to initialize the
clusters. The batch can be small but enough to select $k\cdot p$ medoids. In
the experiments, we used a batch size of ($1.5\cdot k\cdot p$). We use a non-
uniform sampling strategy, similar to the Lloyd’s algorithm [16], to select
the primary medoid of each cluster: SECLEDS selects the first medoid of the
first cluster ($m_{1,1}^{1}$) arbitrarily from the batch. Another $k-1$
medoids are sampled with a probability proportional to the squared distance
between $m_{1,1}^{1}$ and other items in $\mathbb{B}$. This initializes the
primary medoid of each cluster. The other p-1 medoids for each cluster
$C_{i}^{1}$ are the items in $\mathbb{B}$ that are closest to its primary
medoid $m_{i,1}^{1}$. This way, the medoids maintain cluster separation by
reducing the risk of medoids from multiple clusters overlapping each other.
All medoids start with 0 votes.

Assign and Update. With the clusters initialized, the stream processing
begins. Assign and Update are called for each item in $S$. Assign has 3 steps:
(i) An incoming item $s$ at time $t$ is assigned to the cluster $C_{cid}^{t}$
for which its previous medoids $C_{cid}^{t-1}$ have the least average distance
to $s$, formally defined in Eq. (1) for any given distance function $d(.,.)$.
(ii) The closest medoid to $s$ receives a vote, while exponential decay
$\lambda$ is applied to the other medoids _i.e.,_ for all $1\leq j\leq p$ and
$j^{\prime}\neq j$: $v_{cid,j}^{t}=(v_{cid,j}^{t-1}+1)$ if
$j=\operatorname*{arg\,min}_{j}d(s,m_{cid,j}^{t-1})$, otherwise
$v_{cid,j^{\prime}}^{t}=v_{cid,j^{\prime}}^{t-1}\cdot(1-\lambda)$. This way,
the medoids maintain an estimate of their centers of mass without storing any
part of the stream. (iii) The votes of all other clusters remain the same,
_i.e.,_ $v_{i,j}^{t}=v_{i,j}^{t-1}$ for all $i\not=cid$ and $1\leq j\leq p$.

$C_{cid}^{t}=\operatorname*{arg\,min}_{1\leq cid\leq
k}\frac{\sum_{j=1}^{p}d(s,m_{cid,j}^{t-1})}{p}$ (1)

At every time step $t$, the new data item $s$ is promoted to be a medoid of
$C_{cid}^{t}$: the medoid having the least votes which is not the newest
medoid is replaced by $s$, _i.e.,_ $\\{m_{cid,j}^{t}=s,\ v_{cid,j}^{t}=0\\}$
where $j=\operatorname*{arg\,min}_{j}v_{cid,j}^{t-1}$ and
$m_{cid,j}^{t-1}\neq\eta_{cid}^{t}$, where $\eta_{cid}^{t}$ keeps track of the
newest medoid of cluster $cid$ at time $t$. Inspired by Tabu search [11],
including $\eta_{cid}^{t}$ ensures that the most-recently updated medoid is
not selected to be replaced each time. Tabu search is a local search meta-
heuristic that selects which values to change except for the last $\delta$
ones ($\delta=1$ in this case).

Input: Data stream, nclusters, nprototypes: $S,\ k,\ p$

1 function $\textsc{SECLEDS}{(}S,k,p{)}$

2 $b\leftarrow 1.5\cdot k\cdot p$

3 $\mathbb{B}$ $\leftarrow$ Collect b items from $S$

4 $\mathcal{C}$ $\leftarrow$ Init($\mathbb{B}$, $k$, $p$) // Init

5 forall  _$s$ in $S[b:]$_ do

6 $cid\ \leftarrow\ \operatorname*{arg\,min}_{1\leq cid\leq
k}\frac{1}{p}\cdot\sum_{j=1}^{p}d(s,m_{cid,j})$ // Assign

7 $j\ \leftarrow\ \operatorname*{arg\,min}_{j}d(s,m_{cid,j})$ for all $1\leq
j\leq p$

8 $v_{cid,j}\leftarrow(v_{cid,j}+1)$, $v_{cid,j^{\prime}}\leftarrow
v_{cid,j^{\prime}}\cdot(1-\lambda)$ for $j^{\prime}\neq j$

9 $j\ \leftarrow\ \operatorname*{arg\,min}_{j}v_{cid,j}$ where
$m_{cid,j}\neq\eta_{cid}$ for all $1\leq j\leq p$ // Update

10 $m_{cid,j}\leftarrow s,\ v_{cid,j}\leftarrow 0$

11 yield $cid$

12

13

14

15function $\textsc{Init}{(}\mathbb{B},k,p{)}$

16 Choose $m_{1,1}\in\mathbb{B}$ arbitrarily. Let
$C_{1}\leftarrow\\{(m_{1,1},0)\\}$

17 for  _$i\leftarrow 2\dots k$_ do

18 Choose $m_{i,1}\in\mathbb{B}$ with probability $d(m_{i,1},m_{1,1})^{2}$,
$m_{i,1}\neq m_{1,1}$

19 Let $C_{i}\leftarrow\\{(m_{i,1},0)\\}$

20

21 for  _$i\leftarrow 1\dots k$_ do

22 $dist$ $\leftarrow$ $d(b,m_{i,1})$ for all $b\in\mathbb{B}$ and $b\neq
m_{i,1}$

23 Choose $\\{m_{i,2}\dots m_{i,p}\\}$ having smallest values in $dist$

24 Update $C_{i}\leftarrow\\{(m_{i,1},0)\dots(m_{i,p},0)\\}$

25

26 return _$\\{C_{1},\dots,C_{k}\\}$_

27

Algorithm 1 SECLEDS for clustering sequences in evolving streams

##### Time complexity.

Given k clusters, p medoids, b batch size, and n items in the stream, SECLEDS
has a time complexity of $\mathcal{O}(\texttt{n})$: SECLEDS selects the first
medoid at random from the initial batch, and then performs b distance
computations to find the other k-1 primary medoids. The rest of the k(p-1)
medoids are also selected using the same distance information. In total, this
requires $\mathcal{O}(\texttt{kb})$ distance computations. For every $s\in S$,
SECLEDS computes the average distance to each cluster, which requires kp
distance computations. Over an entire run, this gives nkp distance
computations. In the Update module, SECLEDS reallocates medoid votes without
any distance computations, making the runtime negligible. Since k and p are
small user-supplied parameters, the overall runtime complexity is
$\mathcal{O}(\texttt{n})$.

##### Space complexity.

After initialization, SECLEDS only stores the p medoids and their votes for
the k clusters. Since these are (small) user-defined parameters, the space
complexity of SECLEDS is $\mathcal{O}(\texttt{1})$.

##  5 Experimental Setup

##### Datasets.

We use three synthetic and a real dataset containing 2D points and univariate
sequences, see Table 1. The data generation process is given in appendix. The
synthetic datasets are released in the SECLEDS code repository.

Blobs: The blob dataset was created using scikit-learn [22]. The dataset
contains $n=100,000$ two-dimensional points ($d=2$), equally distributed in
$k=10$ classes, with varying standard deviations.

Sine-curve: A sine-curve generator was used to create $k=4$ synthetic
univariate sine curves of length $1,250,000$ each, using varying frequency,
phase and error (see appendix). Each curve is partitioned using a non-
overlapping window of length $w=100$ to obtain the experimental dataset. In
total, $n=50,000$ curves are obtained, equally divided across $k=4$ classes.

Sine-curve-drifted: Incremental concept drift is added to the Sine-curve
dataset by shifting the phase of each curve by a factor of $(drift\cdot
c\\_id)$, where $c\\_id$ is the curve index in the stream, and $drift=0.05$.
Note that adding drift to the frequency of the sine curves produces similar
results.

CTU13-9: CTU13 [10] is an open source dataset composed of network traffic
(netflows) coming from real botnet-infected hosts and normal hosts. We use
scenario-9, containing 10 bot-infected hosts and 6 benign hosts. A total of
$2,087,509$ (normal and botnet) netflows were captured over 5 hours and 37
minutes. We obtain $n=213,386$ univariate sequences of length $w=100$ using a
sliding window model [33] with $step\\_size=1$ (see appendix).

We use Euclidean distance for all datasets. For the Sine-curve and network
traffic datasets, we additionally use Dynamic time warping (DTW). Note that
Euclidean distance can only be used with fixed-length sequences and often
produces less accurate results compared to DTW [29, 12].

Table 1: Summary of experimental datasets.

Dataset | Type | Drift | Stream size (n) | Clusters (k) | Dimensions (d,w)
---|---|---|---|---|---
Blobs | Synthetic | No | 100,000 | 10 | (2,1)
Sine-curve | Synthetic | No | 50,000 | 4 | (1,100)
Sine-curve-drifted | Synthetic | Yes | 50,000 | 4 | (1,100)
CTU13-9 | Real-world | Yes | 213,386 | 2 | (1,100)

##### Stream configuration.

A data stream $S$ of size $n$ is constructed from a chosen experimental
dataset. For each experiment, the clustering task is executed $trials$-times,
randomly shuffling the stream each time, to make the results data-order-
invariant. A clustering task invokes SECLEDS and the baselines such that each
algorithm observes the exact same order of data arrival. In this paper, we set
$trials=10$, unless otherwise reported. All experiments are run on Intel Xeon
E5620 quad-core processor with 74 GB RAM.

##### Evaluation.

We use two metrics for performance evaluation: i) runtime to cluster a stream
size of $n$; ii) F1 score computed from the pairwise co-occurrences of items
in the stream using Eq. (2), as originally defined in [19].

$\displaystyle eval(a,b)=\begin{cases}y_{a}=y_{b}\land C_{x}=C_{y},&\text{true
positive}\\\ y_{a}=y_{b}\land C_{x}\neq C_{y},&\text{false negative}\\\
y_{a}\neq y_{b}\land C_{x}=C_{y},&\text{false positive}\\\ y_{a}\neq
y_{b}\land C_{x}\neq C_{y},&\text{true negative}\\\ \end{cases}$ (2)

where $y_{a}$ and $y_{b}$ are labels of items $a$ and $b$ that are placed in
clusters $C_{x}$ and $C_{y}$. Since clusters do not have pre-defined labels,
data from one class may be assigned to arbitrary clusters in different runs.
Thus, instead of looking at the predicted label, we measure F1 using the
pairwise co-occurrences of true labels.

##### Baselines.

We compare SECLEDS with state-of-the-art open-source partition-based
clustering algorithms with k-fixed clusters: a) Streaming: CluStream,
StreamKM++; b) Batching: MiniBatch k-means; c) Offline: BanditPAM. MiniBatch
k-means and StreamKM++ are online versions of k-means, while CluStream is an
adaptive, online-offline algorithm. BanditPAM (v1.0.5) is used as a benchmark
for the best achievable clustering on a static dataset. We set time_window=1,
max_micro_clusters=$k\cdot p$, halflife=0.5 for CluStream; chunk_size=1,
halflife=0.5 for StreamKM++; batch_size=1, max_iter=1 for Minibatch k-means.

Figure 3: Clustering Blobs and Sine-curves: SECLEDS’s runtime grows
approximately linearly with stream size, while maintaining competitive F1
score with the best-performing baselines, _i.e.,_ BanditPAM and Minibatch
k-means. SECLEDS consistently performs better than all baselines in the
presence of concept drift.

##  6 Empirical Results

Key findings.  In this section, we empirically demonstrate the following
results:

  1. 1.

SECLEDS produces high-quality clusters, regardless of concept drift, stream
size $n$, data dimensionality $(d,w)$, and number of clusters $k$. SECLEDS
shows competitive F1 compared to the best performing baseline (BanditPAM),
while reducing the number of required distance computations by 83.7%.

  2. 2.

SECLEDS outperforms all baselines by 138.7% when the stream contains concept
drift. SECLEDS outperforms the best-performing streaming baseline by 58.2% on
Blobs, 33.3% on Sine-curve, and 143.7% on Sine-curve-drifted.

  3. 3.

SECLEDS-dtw clusters $\sim$5.5h of network traffic in just 8% of the time.
Thus, it can handle networks with bandwidths of up to 1.08 Gbps, which is
significantly higher than the requirements of a typical enterprise network.

##### Point vs. sequence clustering.

We use Blobs with $k=10$ on stream sizes $n=(5000,\dots 50,000)$, and Sine-
curve with $k=3$ on stream sizes $n=(2000,\dots 30,000)$. For both, we set
$p=5$, $\lambda=0.1$, $trials=10$. The mean and standard deviation of the F1
scores and runtimes are given in Figure 3. The benchmark (BanditPAM) achieves
a mean F1 of 0.95 and 1.0 for Blobs and Sine-curve, respectively.

SECLEDS outperforms both CluStream and StreamKM++ on the Blobs dataset, and
additionally outperforms Minibatch k-means on the Sine-curve dataset.
Minibatch k-means performs exceptionally well on point clustering, but loses
its edge on sequence clustering. This is because the centroids are computed by
collapsing temporally-linked dimensions into single values that do not
adequately represent the sequences. An improvement in F1 score is observed for
CluStream and StreamKM++ on the higher dimensional Sine-curve dataset because
of fewer clusters ($k=10$ vs. $k=3$).

We also compare the effect of euclidean and dynamic time warping distance on
the Sine-curve dataset. Although, they both produce equivalent results, it
must be noted that euclidean distance only works with fixed-length sequences.
An example of SECLEDS-dtw on clustering bivariate sequences $d=2,w=$(min:15,
max:121) from UJI Pen Characters [8] is given in the appendix.

##### Initialization quality.

Stream clustering algorithms are greatly impacted by the quality of cluster
initialization. To test this, we compare SECLEDS against SECLEDS-rand
(initialized with randomly selected medoids from the initial batch
$\mathbb{B}$). Evidently, the clusters take a long time to converge,
regardless of the stream size. The cumulative F1 score over time for these
configurations is given in the appendix, showing that although the impact of
poor initialization is reduced over time, SECLEDS-rand does not completely
recover from it. Thus, the distance-based non-uniform sampling strategy proves
to be extremely helpful in initializing good clusters.

##### Clustering with Concept drift.

We use Sine-curve-drifted with $k=3$, $p=5$, $\lambda=0.1$, $trials=10$ on
stream sizes $n=(2000,\dots 30,000)$. SECLEDS outperforms all baselines by
138.7%, and outperforms the best-performing streaming baseline by 143.7%, on
average. BanditPAM no longer serves as a benchmark because it only has a
static view of the data, _i.e.,_ it does not distinguish between class
distributions at $T=t_{x}$ and $T=t_{x+y}$. Both SECLEDS and CluStream
maintain their F1 scores with concept drift, but SECLEDS is 161.8% better than
CluStream. StreamKM++ and Minibatch k-means observe a significant reduction in
their performance. We hypothesize that it might be due to the lack of
exponential decay in k-means, which limits the movement of the centroids
towards newer data. This experiment provides strong evidence for SECLEDS’
ability to handle concept drift with only k-fixed clusters.

##### Runtime analysis.

StreamKM++ and Minibatch k-means are among the fastest clustering algorithms
on all datasets, which is expected since they are based on k-means. CluStream
does not scale well for high-dimensional datasets, and is much slower than
SECLEDS on sequence clustering. As the stream size $n$ grows, SECLEDS also
becomes faster than the high-performance implementation of BanditPAM on both
point and sequence clustering. Interestingly, the runtimes of BanditPAM,
CluStream and StreamKM++ seem to be affected by concept drift: given the same
dataset and constant parameters, their runtimes increase approximately twofold
when there is drift in the data. We hypothesize that this is a side effect of
the sampling strategy used to speed up these algorithms.

(a) $p=\\{3,5\\}$ provide the best trade-off between runtime and F1 for Sine-
curve.

(b) SECLEDS requires 83.7% fewer distance computations compared to BanditPAM.

Figure 4: Scaling with $n$, $k$, and $p$: (a) Empirical results; (b)
Theoretical estimate.

##### Scaling with $n$, $k$, and $p$.

We use Sine-curve with $k=4$, $p=\\{1,3,5,10\\}$, $\lambda=0.1$, $trials=10$
on stream sizes $n=(2,000,\dots 30,000)$. The mean and standard deviation of
the F1 and runtime of SECLEDS is reported in Figure 4(a). A single medoid per
cluster, which is standard for PAM-based algorithms, does poorly in a
streaming setting. Intuitively, more medoids help to improve the stability of
the clusters, but the relationship is not linear. If p is set too low, the
medoids keep jumping to various regions in the dataset, and if it is set too
high, the medoids slow down the evolution of the clusters, having an equally
detrimental effect on the performance. The optimal value of p with respect to
performance and runtime is dataset-dependent. For Sine-curve, $p=\\{3,5\\}$
are good alternatives. Additionally, although SECLEDS has multiple medoids per
cluster, it performs significantly fewer distance computations compared to the
(almost linear) BanditPAM. Figure 4(b) shows this for increasing stream size
$n$, number of clusters $k$, and number of medoids $p$, with BanditPAM as
reference.

Table 2: Clustering real network traffic: Compared to BanditPAM, SECLEDS
requires fewer distance computations, is faster, and has a better cluster
quality. SECLEDS-dtw is slower but produces better clusters than SECLEDS.
Overall, $\sim$5.5h of network traffic is clustered in under 27 minutes (Bold
= best scores).

|  | Stream
---
config.
| # Distances
---
($k=2$)
| Run time
---
($k=2$)
| F1
---
($k=2$)
| F1
---
($k=5$)
BanditPAM |  | Time-ordered
---
$10.3\times 10^{6}$ | 978.03s | 0.64 | 0.38
| Cross-validated
---
984.8s | 0.64 | 0.38
SECLEDS |  | Time-ordered
---
$\mathbf{2.1\times 10^{6}}$ | 629.39s | 0.85 | 0.82
| Cross-validated
---
631.84s | 0.79 | 0.76
SECLEDS-dtw |  | Time-ordered
---
$\mathbf{2.1\times 10^{6}}$ | 1623.05s | 0.85 | 0.88
| Cross-validated
---
1626.89s | 0.81 | 0.80

(a) BanditPAM

(b) SECLEDS-dtw

(c) SECLEDS

Figure 5: Visualizing the medoids of BanditPAM, SECLEDS & SECLEDS-dtw on
$k=2,\ p=5$. Each row is a medoid. The label denotes curve identifier and
$y_{i}$.

###  6.1 Use case: Intelligent network traffic sampling via SECLEDS

A typical enterprise network has a bandwidth of 25
Mbps222https://mosaicnetworx.com/it-challenges/bits-bytes-understanding-
enterprise-network-speeds/, which produces about 17,000 packets per second,
consuming 2 terabytes of storage space each day! To conserve space, the
packets are aggregated into network flows (netflows) at the router level, and
only a fraction of them are stored for analysis _i.e.,_ 1 in $N$ netflows are
stored. Naturally, randomly sampled network traffic does not preserve the
temporal patterns of the data, thus limiting the efficacy of traffic profiling
and behavior analytics.

We propose to cluster sequences of netflows using SECLEDS, and to periodically
store a medoid snapshot of each cluster, since they are representative of the
network traffic seen so far. This way, each snapshot stores an overview of
temporally-linked netflows. The number of clusters $k$ can be chosen depending
on the required granularity of behaviors captured by the clusters. It can also
be approximated from an initial batch using, _e.g.,_[3]. The number of medoids
$p$ can be configured according to available storage space, network bandwidth,
and the intervals at which to store the medoids.

We demonstrate this use case by generating a stream from the CTU13-9 netflows.
The construction and feature engineering processes are given in the appendix.
In short, the ground truth provides two classes, _i.e.,_
$y_{i}\in\\{botnet,normal\\}$. Univariate sequences of average bytes per
netflow are used to separate the two classes. We use two configurations for
the stream: i) Time-ordered: the sequences arrive in order of their
timestamps; ii) Cross-validated: we shuffle the stream. We run SECLEDS and
SECLEDS-dtw with $k=2$, $p=5$, $trials=5$, and compare the performance against
BanditPAM. The results are given in Table 2.

SECLEDS is faster and produces better medoids compared to BanditPAM. Figure 5
visualizes the final medoids produced by all three algorithms in the form of
temporal heatmaps. Temporal heatmaps have previously been used to visualize
temporal similarities in [20]. Each row shows a sequence (medoid), and the
colors indicate the magnitude of the curve at each time step. Both medoids of
BanditPAM are from the normal class. Although the medoids of SECLEDS-dtw are
all from the botnet class, it is evident that they capture distinct behaviors
of the malicious hosts. SECLEDS finds medoids from both classes, but the
clusters are impure, _i.e.,_ more medoids on average are from different
classes. As such, the cluster quality of SECLEDS-dtw is significantly better
than that of SECLEDS.

The results indicate that there are many smaller classes in the network
stream, reflecting the various behaviors of benign and infected hosts. When
$k$ is set to a larger number, the clustering algorithms find smaller, purer
data regions, _e.g.,_ for $k=5$, SECLEDS-dtw produces 4 pure clusters (2
normal and 2 botnet), while SECLEDS only produces 1 pure (normal) cluster, see
appendix for their temporal heatmaps. Table 2 shows the F1 scores for $k=5$.
Note that although the clustering results for $k=5$ are better than $k=2$, the
former obtains a lower F1 score as a side-effect of the metric: it penalizes
higher number of clusters when less class labels are available by lowering the
recall. As such, we recommend to over-estimate $k$ in order to sample many
regions from the network traffic.

Finally, SECLEDS is faster than SECLEDS-dtw, as expected. SECLEDS clusters the
entire stream in 3.1% of the traffic collection time, and SECLEDS-dtw in 8% of
the collection time. This experiment provides evidence that SECLEDS can handle
much larger network bandwidths.

##### Network bandwidth support.

SECLEDS-dtw can handle networks with a bandwidth of up to 1.08 Gbps, which is
more than sufficient for small to medium enterprises. The experiments in Table
2 show that SECLEDS spends $\frac{1626.89}{213386}=0.0076$ seconds on average
to cluster a single sequence of length $w=100$. Thus, SECLEDS can cluster
$131.58$ sequences per second. Assuming that the sequence windows $w$ are non-
overlapping over the traffic stream, SECLEDS can process $13,158$ individual
netflows per second. The CTU13-9 dataset is composed of $115,415,321$ packets
aggregated into $2,087,509$ netflows. Assuming uniform distribution, each
netflow contains about $55.2$ packets. SECLEDS can, thus, process $726,315.79$
packets per second. Given that each network packet is about $1500$ bytes, this
makes a total of $1.089$ Gigabytes per second. Similarly, SECLEDS can handle
network bandwidths of up to 2.79 Gbps.

##  7 Conclusions

We propose SECLEDS, a streaming version of k-medoids with constant memory
footprint. SECLEDS uses a combination of multiple medoids per cluster and a
medoid voting scheme to create k-clusters that evolve with evolving data
streams. Testing on several real and synthetic datasets and comparing against
state-of-the-art baselines, we demonstrate that i) SECLEDS achieves
competitive F1 score compared to the benchmark (BanditPAM) on streams without
concept drift; ii) SECLEDS outperforms all baselines by 138.7% on streams with
concept drift; iii) SECLEDS reduces the number of required distance
computations by 83.7% compared to the benchmark, making it faster than
BanditPAM and CluStream for several clustering tasks, iv) SECLEDS can support
high-bandwidth network streams of up to 1.08 Gbps using the expensive dynamic
time warping distance. These results reinforce the importance of designing
lightweight medoid-based stream clustering algorithms.

####  7.0.1 Acknowledgements

We thank Ruben te Wierik, Silviu Fucarev, and Rami Al-Obaidi for their
contributions to the SECLEDS algorithm.

## References

  * [1] Ackermann, M.R., Märtens, M., Raupach, C., Swierkot, K., Lammersen, C., Sohler, C.: Streamkm++ a clustering algorithm for data streams. JEA 17, 2–1 (2012)
  * [2] Aggarwal, C.C., Philip, S.Y., Han, J., Wang, J.: A framework for clustering evolving data streams. In: VLDB. pp. 81–92. Elsevier (2003)
  * [3] de Andrade Silva, J., Hruschka, E.R.: Extending k-means-based algorithms for evolving data streams with variable number of clusters. In: ICMLA. vol. 2, pp. 14–19. IEEE (2011)
  * [4] Barros, R.S.M., Santos, S.G.T.C.: A large-scale comparison of concept drift detectors. Information Sciences 451, 348–370 (2018)
  * [5] Boeva, V., Nordahl, C.: Modeling evolving user behavior via sequential clustering. In: ECML-PKDD. pp. 12–20. Springer (2019)
  * [6] Cao, F., Estert, M., Qian, W., Zhou, A.: Density-based clustering over an evolving data stream with noise. In: SDM. pp. 328–339. SIAM (2006)
  * [7] Cook, D.J., Krishnan, N.C., Rashidi, P.: Activity discovery and activity recognition: A new partnership. IEEE transactions on cybernetics 43(3), 820–828 (2013)
  * [8] Dua, D., Graff, C.: UCI machine learning repository (2017)
  * [9] Fahy, C., Yang, S.: Finding and tracking multi-density clusters in online dynamic data streams. IEEE Transactions on Big Data (2019)
  * [10] Garcia, S., Grill, M., Stiborek, J., Zunino, A.: An empirical comparison of botnet detection methods. computers & security 45, 100–123 (2014)
  * [11] Glover, F.: Future paths for integer programming and links to artificial intelligence. Computers & operations research 13(5), 533–549 (1986)
  * [12] Guijo-Rubio, D., Durán-Rosal, A.M., Gutiérrez, P.A., Troncoso, A., Hervás-Martínez, C.: Time-series clustering based on the characterization of segment typologies. IEEE transactions on cybernetics 51(11), 5409–5422 (2020)
  * [13] Guo, J., Liu, G., Zuo, Y., Wu, J.: Learning sequential behavior representations for fraud detection. In: ICDM. pp. 127–136. IEEE (2018)
  * [14] Hyde, R., Angelov, P., MacKenzie, A.R.: Fully online clustering of evolving data streams into arbitrarily shaped clusters. Information Sciences 382, 96–114 (2017)
  * [15] Islam, M.K., Ahmed, M.M., Zamli, K.Z.: A buffer-based online clustering for evolving data stream. Information Sciences 489, 113–135 (2019)
  * [16] Lloyd, S.: Least squares quantization in pcm. IEEE transactions on information theory 28(2), 129–137 (1982)
  * [17] Lu, J., Liu, A., Dong, F., Gu, F., Gama, J., Zhang, G.: Learning under concept drift: A review. TKDE 31(12), 2346–2363 (2018)
  * [18] Lu, N., Zhang, G., Lu, J.: Concept drift detection via competence models. Artificial Intelligence 209, 11–28 (2014)
  * [19] Manning, C., Raghavan, P., Schütze, H.: Introduction to information retrieval. Natural Language Engineering 16(1), 100–103 (2010)
  * [20] Nadeem, A., Hammerschmidt, C., Gañán, C.H., Verwer, S.: Beyond labeling: Using clustering to build network behavioral profiles of malware families. In: MAAIDL, pp. 381–409. Springer (2021)
  * [21] Nadeem, A., Verwer, S., Moskal, S., Yang, S.J.: Alert-driven attack graph generation using s-pdfa. IEEE Transactions on Dependable and Secure Computing 19(2), 731–746 (2021)
  * [22] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., et al.: Scikit-learn: Machine learning in python. Journal of machine Learning research 12, 2825–2830 (2011)
  * [23] Schubert, E., Rousseeuw, P.J.: Faster k-medoids clustering: improving the pam, clara, and clarans algorithms. In: SISAP. pp. 171–187. Springer (2019)
  * [24] Sculley, D.: Web-scale k-means clustering. In: WWW. pp. 1177–1178 (2010)
  * [25] Silva, J.A., Faria, E.R., Barros, R.C., Hruschka, E.R., Carvalho, A.C.d., Gama, J.: Data stream clustering: A survey. CSUR 46(1), 1–31 (2013)
  * [26] Tiwari, M., Zhang, M.J., Mayclin, J., Thrun, S., Piech, C., Shomorony, I.: Banditpam: Almost linear time k-medoids clustering via multi-armed bandits. NeurIPS 33, 10211–10222 (2020)
  * [27] Ushakov, A.V., Vasilyev, I.: Near-optimal large-scale k-medoids clustering. Information Sciences 545, 344–362 (2021)
  * [28] Wang, T., Li, Q., Bucci, D.J., Liang, Y., Chen, B., Varshney, P.K.: K-medoids clustering of data sequences with composite distributions. IEEE Transactions on Signal Processing 67(8), 2093–2106 (2019)
  * [29] Wang, X., Mueen, A., Ding, H., Trajcevski, G., Scheuermann, P., Keogh, E.: Experimental comparison of representation methods and distance measures for time series data. Data Mining and Knowledge Discovery 26(2), 275–309 (2013)
  * [30] Wang, Y., Chen, L., Mei, J.P.: Incremental fuzzy clustering with multiple medoids for large data. IEEE transactions on fuzzy systems 22(6), 1557–1568 (2014)
  * [31] Zhang, T., Ramakrishnan, R., Livny, M.: Birch: A new data clustering algorithm and its applications. Data mining and knowledge discovery 1(2), 141–182 (1997)
  * [32] Žliobaitė, I., Pechenizkiy, M., Gama, J.: An overview of concept drift applications. Big data analysis: new algorithms for a new society pp. 91–114 (2016)
  * [33] Zubaroğlu, A., Atalay, V.: Data stream clustering: a review. Artificial Intelligence Review 54(2), 1201–1236 (2021)

##  Appendix 0.A Dataset generation process

Figure 6: Experimental datasets: (Top left): Scatter plot showing the Blobs
dataset ($k=10$). (Top middle): t-SNE plot showing a 2D representation of the
Sine-curve dataset ($k=4$). (Top right): t-SNE plot showing a 2D
representation of the CTU13-9 dataset ($k=2$). (Bottom left): Temporal heatmap
showing 20 randomly sampled sequences from Sine-curve. (Bottom middle)
Temporal heatmap showing 20 sequences from Sine-curve-drifted. (Bottom right):
Temporal heatmap showing 10 sequences from CTU13-9.

Sine-curve.  The parameters used to create the 4 sine curves are:
$frequency=\\{(0.1,0.12),(0.2,0.22),(0.4,0.42),(0.6,0.62)\\}$,
$error=\\{0.2,0.4,0.7,0.1\\}$ and $phase=\\{5,12,-10,-20\\}$. For each curve
$i\in\\{1,\dots k\\}$: (a) the phase value is $phase_{i}$, (b) a random
frequency is selected from $frequency_{i}$ range, and (c) a random error of
factor $error_{i}$ is added at every time step of the curve. Figure 6 shows
the 2 dimensional representation of the dataset via a t-SNE plot and temporal
heatmap. 5 randomly sampled sequences from the 4 classes are shown in the
heatmap. The x-axis represents time, thus each row is a sequence. The colors
show the magnitude of the curve at each time step. Temporal heatmaps have
previously been used to visualize temporal similarities in (Nadeem,2021).

CTU13-9.  The network traffic is divided into several sub-streams based on the
source IP address. There are 16 sub-streams, each collecting netflows for one
of the 16 hosts in the dataset (10 infected with Neris botnet, 6 benign). We
assume that the netflows arrive in sequences of length $w=100$, which are
constructed by sliding a window $w$ over each of the sub-streams with
$stepsize=1$. This results in $n=213,386$ netflow sequences. Using the two
natural classes in the dataset, we follow (Garcia,2014) and label the
sequences coming from infected sub-streams as botnet, while the others as
normal.

A preliminary feature engineering is conducted to find the netflow feature(s)
that differentiates the two classes. There is extensive research studying the
impact of various features in traffic classification. For simplicity, we
consider a single feature at a time. Searching for better feature combinations
is left as future work. We visualize the class separation obtained by i)
netflow duration, ii) total bytes, iii) total packets, iv) inter-arrival time
between netflows, v) average bytes per netflow, and vi) average packets per
netflow. We select average bytes to separate the two classes. Figure 6 shows a
t-SNE plot of the sequences. The unclear class separation emphasizes just how
challenging the clustering task is.

##  Appendix 0.B Cluster initialization quality

(a) Blobs, $k=10;p=5;n=5,000$

(b) Blobs, $k=10;p=5;n=20,000$

Figure 7: F1 score over time for SECLEDS and SECLEDS-rand. Randomly
initialized clusters perform sub-optimally regardless of stream size.

Figure 7 shows the F1 scores over time for SECLEDS and SECLEDS-rand (randomly
initialized medoids) for two stream sizes. It is evident that poorly
initialized clusters are unable to perform optimally. Since the medoids in a
single cluster are blind to each other, once poorly initialized, the clusters
never converge. Further investigation is required to make the medoids aware of
each other since the additional distance computations will increase the
complexity of SECLEDS.

##  Appendix 0.C Network traffic sampling with SECLEDS

The CTU13-9 dataset contains two big classes corresponding to the normal and
botnet netflows. However, the experiments with $k=2$ indicate that there are
several smaller regions in the data distribution, reflecting to the various
behaviors of the normal and bot-infected hosts. We demonstrate this by setting
$k=5$. Similar to the $k=2$ case, BanditPAM finds all clusters in the normal
class. Figure 8 shows the final medoids obtained by SECLEDS-dtw and SECLEDS.
It shows that SECLEDS-dtw finds four fully pure clusters, _i.e.,_ two in the
botnet class and two in the normal class. SECLEDS (with Euclidean distance)
struggles to find more than one pure cluster. This experiment provides
evidence for alignment-based distance measures, such as dynamic time warping,
to have superior performance than Euclidean distance for sequence clustering
tasks.

(a) SECLEDS-dtw

(b) SECLEDS

Figure 8: Visualizing the medoids of SECLEDS-dtw & SECLEDS on $k=5,\ p=5$.
Each row is a medoid. The label denotes curve identifier and $y_{i}$. The pure
clusters are highlighted by green for normal and red for botnet class.

##  Appendix 0.D Clustering Handwritten character strokes

Despite the growing interest in sequential features over the years, one of the
biggest bottlenecks in clustering sequences is the choice of an apt distance
measure. Euclidean distance is by far the most popular metric used for
clustering tasks. However, Euclidean distance does not perform accurately for
out-of-sync sequential data, which is a common property of sequences. Although
expensive, alignment-based distances, such as dynamic time warping, have
instead been shown to be more accurate. Another benefit of dynamic time
warping is that it can handle sequences of arbitrary lengths. We demonstrate
this by clustering a dataset of handwritten characters from the UCI machine
learning repository (Dua,2017). The UJI Pen Characters dataset contains 467
bivariate sequences from 13 classes, _i.e.,_
$\\{C,O,S,U,V,W,1,2,3,5,6,8,9\\}$. The length of the shortest character stroke
is 15, while that of the longest stroke is 121. Figure 9 shows the t-SNE plot
of the dataset.

SECLEDS-dtw is the only algorithm that can cluster this dataset straight out-
of-the-box. The other baselines use Euclidean distance, which does not support
arbitrary-length sequences. Figure 9 also shows the 5-fold cross-validated F1
score over time for clustering the character strokes. BanditPAM is used as a
reference benchmark, which requires the pairwise distance matrix of all
sequences to cluster the dataset. Note that such a distance matrix is not
available in a streaming setting, since it requires the entire dataset to be
stored in memory. Additionally, we compare the performance with a randomly
initialized SECLEDS, referred to as SECLEDS-dtw-rand. BanditPAM obtains the
highest achievable F1 score of 0.294 since very few classes are properly
separated. SECLEDS-dtw obtains a comparable F1 score of 0.291, while SECLEDS-
dtw-rand obtains a score of 0.247. This experiment shows that dynamic time
warping is indeed a more suitable distance measure for sequence clustering.

Figure 9: Clustering the UJI Pen Characters dataset: (Left): t-SNE plot
showing the 2D representation of the character strokes. (Right): F1 score over
time for BanditPAM, SECLEDS-dtw and SECLEDS-dtw-rand. SECLEDS-dtw has
comparable F1 score to BanditPAM, demonstrating that dynamic time warping is
more suitable for sequence clustering in a streaming setting.


# Orbit-localised thermosphere density prediction using a Kalman filter based
calibration of empirical models

George Bowden

###### Abstract

Accurate estimation of thermosphere mass density is critical to determining
how satellite orbits evolve over time and thus to planning and managing space
missions. Empirical thermosphere models are commonly employed for this
purpose, but have substantial uncertainties. In this work, a Kalman filter
method for calibrating these models along a particular satellite’s trajectory
is described. This method was applied to calibrate the NRLMSISE-00, JB2008,
and DTM-2020 models with respect to densities measured by the Swarm-C
satellite. Substantial improvements in root mean squared density residuals
were obtained using the technique when compared with either uncalibrated model
output or calibration using a linear regression on previous data. Further
improvement was obtained by combining estimates from different models using a
best linear unbiased estimator method.

###### keywords:

Thermosphere density , Modeling , Kalman filters

###### PACS:

92.60.hb , 95.40.+s

###### MSC:

93E11

††journal: Acta Astronautica

[inst1]organization=School of Engineering and Information
Technology,addressline=The University of New South Wales, city=Canberra,
postcode=2612, state=ACT, country=Australia

A method is described for estimating the time-varying linear calibration of
thermosphere density models for a particular satellite using a Kalman
filtering algorithm.

The method is applied to orbit means and principal components of densities
from precise orbit determination of the Swarm C satellite.

##  1 Introduction

Aerodynamic drag is the greatest source of uncertainty for the determination
and propagation of satellite trajectories in Low Earth Orbit (LEO). This force
arises from interaction with the thermosphere, the atmospheric layer from
$\sim 90$ to $600$ km altitude. Accurately estimating satellite drag is a key
challenge for mission planning, lifetime estimation, reentry prediction,
collision avoidance, and propulsionless manoeuvring. Satellite drag, $F_{D}$,
can be expressed in terms of atmospheric density, $\rho$, drag coefficient,
$C_{D}$, speed relative to the atmosphere, $v$, and reference area, $A$,

$F_{D}=\frac{C_{D}\rho v^{2}A}{2}.$ (1)

For a particular satellite, $F_{D}$ may be inferred from on-board
accelerometer measurements or orbit determination from some combination of
radar, laser ranging, and Global Navigation Satellite System (GNSS)
measurements.

Significant uncertainty typically exists regarding thermospheric density and
satellite aerodynamic properties [1]. Drag coefficient is determined by
complex flow physics and therefore challenging to compute accurately. It is
dependent on often complicated satellite geometries, gas-surface interactions,
and thermosphere composition, though is not sensitive to density, velocity,
and length scale. Thermosphere densities vary due to fluctuations in incoming
extreme ultraviolet (EUV) radiation, electrical energy and energetic particles
from the solar wind, and waves and instabilities propagating from the lower
atmosphere [2]. Empirical thermosphere models incorporate the former two
influences through parameterisations based on climatological data. However,
their accuracy is limited by the relative sparseness of underlying data and
unmodelled physics.

An alternative to applying global thermosphere models to individual satellites
is to obtain orbit-localised density estimates for those satellites based on
their previous density measurements. Here density estimates are restricted to
individual satellite trajectories and based on prior density measurements from
the satellites. An example of such an approach, employing a linear
autoregressive filter to forecast density data from the CHAllenging
Minisatellite Payload (CHAMP) satellite, was described by Stastny et al. [3].
Similarly, Pérez et al. detailed the use of feed-forward time-delay artificial
neural networks (ANNs) to predict density for CHAMP [4]. Solar and geomagnetic
index data were included as inputs to incorporate space weather effects. These
computationally inexpensive models could feasibly be implemented as part of
on-board orbit propagation.

In addition to the aforementioned work, Pérez and Bevilacqua investigated the
use of ANNs for calibrating density predictions from empirical thermosphere
models [5]. Their networks were trained and validated on the Gravity Recovery
and Climate Experiment (GRACE) A data, though tested on CHAMP data. The
influence of uncertainty in the drag coefficient of a satellite on its orbit
determination and propagation could be reduced by calibrating density
estimates for that particular satellite. This calibration would help ensure
that density and drag coefficient estimates together estimate drag
accelerations accurately, despite there being significant uncertainty in the
former quantities individually.

Here a Kalman filter based method for calibrating thermosphere density
estimates from empirical models along the trajectory of a particular satellite
is described and evaluated. Such a method can be used with irregularly spaced
and delayed density data, by contrast with prior orbit-localised methods. The
method can account for data with non-uniform uncertainties. It also has the
advantage of providing uncertainties associated with its estimates.
Furthermore, the filter parameters indicate how rapidly the relationship
between empirical model estimates and satellite drag changes.

In this Kalman filter implementation, the calibration is assumed to be linear
with temporal evolution described by a stochastic process. It is iteratively
updated based on measured orbit densities. This technique is applied to
historical data for the Swarm-C satellite, one component of a constellation
studying Earth’s magnetic field [6].

##  2 Empirical thermosphere models

Empirical thermosphere models represent the dependence of density on
parameters such as altitude, local time, latitude, day of year, solar EUV
flux, and geomagnetic activity. They are considered climatological, being
based on average observed behaviour, though usually incorporate certain
physics-based constraints. The development of such models is discussed in
review articles by Vallado and Finkleman [1] and Emmert [2]. Estimates vary
between models due to use of different thermosphere data-sets in their
derivation, different space weather indices as input, and different fitting
functions to represent dependence on time and location. Three empirical
models, summarised below, are considered in this study. The model
implementations used here are available online [7, 8, 9].

###  2.1 NRLMISE-00

The Naval Research Laboratory Mass Spectrometer Incoherent Scatter radar
Extended 2000 model (NRLMSISE-00) [10] belongs to the MSIS family of empirical
thermosphere models [11, 12]. This model provides mass density, number
densities of major thermosphere species, and neutral temperature. It was
derived from a combination of satellite mass spectrometer (MS), incoherent
scatter radar (ISR), satellite accelerometer, orbit determination,
occultation, and sounding rocket data. Solar EUV variation is included through
daily and $81$-day averaged $F_{10.7}$ solar flux indices (ground-based radio
proxies for EUV radiation) and geomagnetic activity is included through a
history of $3$-hourly planetary $a_{p}$ indices.

###  2.2 JB2008

The Jacchia-Bowman 2008 model (JB2008) is based on satellite orbit
determination and accelerometer data [13]. The model computes mass density and
temperature. As well as the $F_{10.7}$ index, the model incorporates
$S_{10.7}$, $M_{10.7}$, and $Y_{10.7}$ indices derived from satellite-based
EUV, ultraviolet, and X-ray measurements respectively. Geomagnetic activity is
represented through estimated changes in exosphere temperature, which are
derived from the Disturbance Storm Time ($Dst$) index.

###  2.3 DTM-2020

The Drag Temperature Model 2020 (DTM2020) [14] was developed from the DTM
family of empirical thermosphere models [15, 16]. These models are fitted to
satellite MS, ISR, satellite accelerometer, orbit determination, and optical
spectrometer measurements. The version of the model in this study takes as
input daily and $81$-day averaged $F_{10.7}$ to represent EUV variation and
$3$-hour delayed and $24$-hour averaged $K_{p}$ indices to represent
geomagnetic activity.

##  3 Kalman filter approach

The Kalman filter is a widely-used algorithm for recursively estimating the
state of a system based on noisy measurements [17]. This filter determines the
most likely state of a system with additive Gaussian measurement and process
noise. It is applicable to smoothing, filtering, and prediction problems.
Here, the Kalman filter is used to estimate coefficients for calibration of
empirical model orbit mean thermosphere mass density encountered by a
satellite. Scores for each thermosphere mass density principal component can
be calibrated using a similar method. Combining these data allow estimation of
thermosphere mass density at specific locations on the satellite orbit.

###  3.1 Transition model and measurements

Assume that the mean thermosphere mass density over an orbit has an
expectation value which is linearly related to that predicted by an empirical
model, $\rho_{m}$ (i.e. $\left<\rho\right>=m\rho_{m}+c$ where $m$ and $c$ are
constants). Furthermore, assume that the measured value is perturbed by a
stochastic component with a Gaussian distribution
$p\left(\epsilon_{\rho}\right)\sim\mathcal{N}\left(0,R\right)$. This term
represents both measurement noise and the effect of transient perturbations
excluded from climatological models, such as atmospheric waves. Thus, the
measured mean density can be expressed as

$\rho=m\rho_{m}+c+\epsilon_{\rho}.$ (2)

A state vector can be constructed from the coefficients in Eq. 2 for time-step
$k$,

$\mathbf{x}_{k}=\left[\begin{array}[]{c}m_{k}\\\ c_{k}\end{array}\right].$ (3)

Temporal evolution is modelled as a Gaussian random walk process,

$\mathbf{x}_{k+1}=\mathbf{F}_{k}\mathbf{x}_{k}+\mathbf{w}_{k}$ (4)

where state transition matrix $\mathbf{F}_{k}$ is the identity matrix and
$p\left(\mathbf{w}_{k}\right)\sim\mathcal{N}\left(0,\mathbf{Q}_{k}\right)$
where $\mathbf{Q}_{k}=\left(t_{k+1}-t_{k}\right)\mathbf{M}$. Eq. 2 can be
rewritten in matrix form to represent an observation at time $k$,

$\mathbf{z}_{k}=\mathbf{H}_{k}\mathbf{x}_{k}+\mathbf{v}_{k},$ (5)

where $\mathbf{z}_{k}=\rho$, $\mathbf{H}_{k}=\left[\begin{array}[]{c
c}\rho_{m}&1\end{array}\right]$, and
$p\left(\mathbf{v}_{k}\right)\sim\mathcal{N}\left(0,R\right)$.

###  3.2 Filter algorithm

A Kalman filter provides the optimal (maximum likelihood) estimate of the
updated state, $\hat{\mathbf{x}}_{k}$, and its associated covariance matrix,
$\mathbf{P}_{k}$, based on prior values, $\hat{\mathbf{x}}_{k}^{-}$ and
$\mathbf{P}_{k}^{-}$, and observation $\mathbf{z}_{k}$. First, the state and
covariance are updated based on the density measurement at time $k$,

$\hat{\mathbf{x}}_{k}=\hat{\mathbf{x}}_{k}^{-}+\mathbf{K}_{k}\left(\mathbf{z}_{k}-\mathbf{H}_{k}\hat{\mathbf{x}}_{k}^{-}\right)$
(6)
$\mathbf{P}_{k}=\left(\mathbf{I}-\mathbf{K}_{k}\mathbf{H}_{k}\right)\mathbf{P}_{k}^{-},$
(7)

where $\mathbf{K}_{k}$ is the optimal Kalman gain,

$\mathbf{K}_{k}=\mathbf{P}_{k}^{-}\mathbf{H}_{k}^{T}\left(\mathbf{H}_{k}\mathbf{P}_{k}^{-}\mathbf{H}_{k}^{T}+\mathbf{R}_{k}\right)^{-1}.$
(8)

The state is then advanced in time

$\hat{\mathbf{x}}_{k+1}^{-}=\mathbf{F}_{k}\hat{\mathbf{x}}_{k}$ (9)
$\mathbf{P}_{k+1}^{-}=\mathbf{F}_{k}\mathbf{P}_{k}\mathbf{F}_{k}^{T}+\mathbf{Q}_{k}.$
(10)

This provides prior values for the next measurement update, so the process can
be repeated.

In practice, for nowcasting and forecasting such a calibration must be used
some time period after the data from which it was derived. If this is done $n$
time steps ahead, the estimated density and measurement uncertainty are found
to be

$\tilde{\rho}_{k+n}=\mathbf{H}_{k+n}\left(\prod_{i=1}^{n}\mathbf{F}_{k+i-1}\right)\hat{\mathbf{x}}_{k}$
(11)

and

$\tilde{\sigma}_{k+n}=\mathbf{H}_{k+n}\left(\left(\prod_{i=1}^{n}\mathbf{F}_{k+i-1}\right)\mathbf{P}_{k}\left(\prod_{i=1}^{n}\mathbf{F}_{k+i-1}\right)^{T}\right.\\\
+\left.\sum_{i=1}^{n}\left(\prod_{j=i+1}^{n}\mathbf{F}_{k+j-1}\right)\mathbf{Q}_{k+i-1}\left(\prod_{j=i+1}^{n}\mathbf{F}_{k+j-1}\right)^{T}\right)\mathbf{H}_{k+n}^{T}+\mathbf{R}_{k+n}.$
(12)

For the system described in Subsection 3.1 these equations reduce to

$\tilde{\rho}_{k+n}=\mathbf{H}_{k+n}\hat{\mathbf{x}}_{k}$ (13)

and

$\tilde{\sigma}_{k+n}=\mathbf{H}_{k+n}\left(\mathbf{P}_{k}+n\mathbf{M}\right)\mathbf{H}_{k+n}^{T}+R.$
(14)

###  3.3 Parameter estimation

The filter and process noise parameters, $R$ and $\mathbf{M}$, are challenging
to estimate on theoretical grounds. Both short time-scale thermosphere density
perturbations and random errors in density measurements are expected to
contribute significantly to $R$. These may not be independent, given that both
can be associated with changes in thermosphere composition. Moreover, the
somewhat abstract nature of a state vector formed from calibration
coefficients makes it difficult to relate $\mathbf{M}$ to known quantities.
Consequently, the best approach is to determine $R$ and $\mathbf{M}$ based on
filter performance for a set of data. These are referred to as ‘training’
data, by analogy with machine learning.

A maximum likelihood method is used to compute $R$ and $\mathbf{M}$. The
marginal probability density $p\left(R,\mathbf{M}|\rho_{i}\right)$ is
maximised, where $\rho_{i}$ is the set of $m$ density measurements in the
training data. Following the derivation of Eq. 14 of Mehra’s paper on adaptive
filtering [18], the likelihood estimate can be obtained by maximising

$L\left(R,\mathbf{M}\right)=-\frac{1}{2}\sum_{i=1}^{m}\left\\{\frac{\left(\rho_{i}-\mathbf{H}_{i}\hat{\mathbf{x}}_{i-n}\right)^{2}}{\mathbf{H}_{i}\left(\mathbf{P}_{i-n}+n\mathbf{M}\right)\mathbf{H}_{i}^{T}+R}\right.\\\
\left.+\ln\left(\mathbf{H}_{i}\left(\mathbf{P}_{i-n}+n\mathbf{M}\right)\mathbf{H}_{i}^{T}+R\right)\right\\}.$
(15)

Here it is assumed that there is no prior information regarding the values of
$R$ and $\mathbf{M}$, which would otherwise result in the addition of the term
$\ln{\left(p\left(R,\mathbf{M}\right)\right)}$ to the right hand side of the
above equation.

Unconstrained numerical optimisation is performed using the log-Cholesky
parametrisation $M=LL^{T}$ where $L$ is a real lower triangular matrix [19].
The optimisation parameters are the logarithms of on-diagonal elements of $L$
and values of off-diagonal elements. Doing so ensures that $M$ is a positive
definite symmetric matrix, which are essential properties of a realistic
covariance matrix. Similarly, using the optimisation parameter $\ln{R}$ to
represent $R$ ensures that this parameter has a realistic positive value.

###  3.4 Combining density estimates

Using the above, Kalman filters can be obtained for calibrating various
empirical thermosphere models. This leads to the question of how to best
combine their density estimates. Residuals in these predictions cannot be
considered independent, given the common contributions of measurement errors
and non-climatological perturbations. The best linear unbiased estimate based
on correlated estimates is [20]

$\tilde{\rho}_{c}=\boldsymbol{\alpha}^{T}\tilde{\boldsymbol{\rho}},$ (16)

where the weight vector is

$\boldsymbol{\alpha}=\frac{\mathbf{K}_{\boldsymbol{\epsilon}_{\rho}\boldsymbol{\epsilon}_{\rho}}^{-1}\mathbf{u}}{\mathbf{u}^{T}\mathbf{K}_{\boldsymbol{\epsilon}_{\rho}\boldsymbol{\epsilon}_{\rho}}^{-1}\mathbf{u}}$
(17)

with $\mathbf{u}$ being a vector in which each element is unity. The variance
is estimated by

$\sigma_{\rho_{c}}^{2}=\boldsymbol{\alpha}^{T}\mathbf{K}_{\boldsymbol{\epsilon}_{\rho}\boldsymbol{\epsilon}_{\rho}}\boldsymbol{\alpha}.$
(18)

In these equations, $\tilde{\boldsymbol{\rho}}$ is a vector of calibrated
orbit mean densities and
$\mathbf{K}_{\boldsymbol{\epsilon}_{\rho}\boldsymbol{\epsilon}_{\rho}}$ is the
covariance matrix for residuals $\epsilon_{\rho}$.

While Eq. 14 provides an estimate of the variance for the estimates derived
from each empirical model, it does not provide a straightforward method of
computing off-diagonal terms of
$\mathbf{K}_{\boldsymbol{\epsilon}_{\rho}\boldsymbol{\epsilon}_{\rho}}$.
Therefore,
$\mathbf{K}_{\boldsymbol{\epsilon}_{\rho}\boldsymbol{\epsilon}_{\rho}}$ is
instead estimated based on residuals for the training data. As
$\mathbf{P}_{k}$ is found to stabilise very quickly for constant
$\mathbf{Q}_{k}$ and $\mathbf{R}_{k}$ in practice, the assumption of constant
$\mathbf{K}_{\boldsymbol{\epsilon}_{\rho}\boldsymbol{\epsilon}_{\rho}}$ is
reasonable.

##  4 Example case: Swarm-C densities

###  4.1 Swarm-C precise orbit determination densities

The thermosphere mass density data used in this study were values published by
ESA derived from Swarm-C GNSS measurements using Precise Orbit Determination
(POD) [21]. This entailed estimation of aerodynamic accelerations using an
extended Kalman filter method with a force model incorporating gravitational
and radiation forces. Densities were calculated from these accelerations using
direct simulation Monte Carlo simulation to compute $C_{D}$. Thermosphere
composition was taken from NRLMSISE-00 and a constant accommodation
coefficient (ratio of energy exchanged in gas-surface interaction to its
thermodynamic maximum) assumed. The fastest density variations this method
captured were a factor of $4$-$5$ times shorter than the orbital period,
corresponding to a resolution of about $80^{\circ}$ in argument of latitude,
$u$. However, the published data had a $30$ s sampling rate.

###  4.2 Data selection

Density measurements from 2018 and 2019 constituted the training and test data
sets respectively. These data coincided with the solar minimum at the end of
Solar Cycle 24. During this time period, Swarm-C had an altitude of
approximately $450$ km and orbital inclination of $87.4^{\circ}$. Slow nodal
precession occurred such that local time of the ascending node (LTAN) covered
the full $24$ hour range over a period of approximately $9$ months. Training
and test orbit mean density data are plotted along with space weather indices
in Fig. 1.

Figure 1: Swarm-C orbit mean density, $\rho$, geomagnetic activity, $a_{p}$,
and solar flux index, $F_{10.7}$, for $2018$ and $2019$.

###  4.3 Data processing

For each orbit during the training and test periods, Swarm-C POD density data
were interpolated to a grid with spacing $1^{\circ}$ in satellite argument of
latitude, $u$. Orbits with missing data or erroneous values were excluded from
the analysis. Mean orbit densities were then computed averaging over $u$.

Principal component coefficients were computed for the training period data.
This was done without the usual practice of subtracting mean values for each
$u$, as means for a training period will not be applicable for different
periods of the solar cycle. Consequently, the first principal component
reflects these mean values in addition to changes driven by geomagnetic and
solar activity. The first $6$ principal components, those accounting for more
than $1\%$ variation in density, are plotted in Fig. 2. These were used to
compute scores, $t_{i}$, for both the training and test data sets, plotted in
Fig. 3.

The first principal component largely represents the global mean density, with
changes driven by variations in geomagnetic and solar activity. Peaks in its
value correspond to times of elevated $a_{p}$ index in Fig. 2. By contrast,
the second principal component correlates with density changes due to nodal
precession. Fig. 3 indicates that the corresponding score varies with mean $0$
and period of about $9$ months. This function is continuous from
$u=360^{\circ}$ to $0^{\circ}$ and is approximately odd symmetric about
$u=90^{\circ}$, as expected for a circular polar orbit. The third component
seems to reflect seasonal variation in the thermosphere itself. It has period
$12$ months, is continuous from $u=360^{\circ}$ to $0^{\circ}$ and is
approximately odd symmetric about $u=0^{\circ}$. The fifth and sixth principal
components have large discontinuities from $u=360^{\circ}$ to $0^{\circ}$
corresponding to changes in density during the time the satellite takes to
complete an orbit.

Figure 2: Principal components of Swarm-C precise orbit determination density
as functions of argument of latitude, $u$, for $2018$. Coefficients for the
first (blue), second (orange), third (yellow), fourth (purple), fifth (green),
and sixth (cyan) components are plotted. Color available in online version.
Figure 3: Scores for Swarm-C precise orbit determination density in $2018$ and
$2019$, for the principal components plotted in Fig. 2. Scores for first
(blue), second (orange), third (yellow), fourth (purple), fifth (green), and
sixth (cyan) components are plotted.

###  4.4 Orbit mean data

The Kalman filter based calibration technique detailed in Section 3 was
applied to the orbit mean density data and estimates of the NRLMSISE-00,
JB2008, and DTM-2020 models. For the purpose of this example, the NRLMSISE-00
and DTM-2020 empirical models were driven with the definitive solar flux and
geomagnetic indices available from CLS (Collecte Localisation Satellites) [22]
and GFZ (German Research Centre for Geosciences) [23] respectively. The JB2008
model used data available from [24].

The method was applied estimating the calibration coefficient $1$ and $3$ days
ahead. This could represent the case of very short-term density prediction
(nowcasting) where accurate and timely solar flux and geomagnetic indices are
available but satellite-based density measurments are delayed. Alternatively,
it can be considered a best case scenario for density forecasts where solar
flux and geomagnetic indices can be predicted very accurately.

For comparison purposes, fixed calibration coefficients were obtained from
linear least squares regression. These were computed separately from both the
training and the test data. The former represented calibration of test data
based solely on training data. The latter represented the best possible
calibration with fixed coefficients for the test data. However, linear
regression is clearly not a practical method of obtaining this calibration as
it incorporates future data.

The root mean square (RMS) residuals for orbit mean densities estimated using
the four calibration methods are tabulated in Tab. 1. These may be compared
with the mean value of $0.1496\times 10^{-12}$ kg.m-3 for the test data. The
calibrated RMS residuals in the table can be compared with the corresponding
values without calibration; for NRLMSISE-00, JB2008, and DTM-2020 these are
$0.1720\times 10^{-12}$ kg.m-3, $0.0548\times 10^{-12}$ kg.m-3, and
$0.0349\times 10^{-12}$ kg.m-3 respectively. Tab. 1 also includes residuals
for results combined using the method detailed in Subsection 3.4.

Table 1: RMS residuals for different orbit mean density calibration methods. All quantities are shown in $10^{-12}$ kg.m-3 (without brackets) and ratio with the mean value (within brackets). Calibration method | NRLMSISE-00 | JB2008 | DTM-2020 | Combined
---|---|---|---|---
Linear regression |  |  |  |
on training data | 0.0469 (0.31) | 0.0379 (0.25) | 0.0478 (0.32) | 0.0411 (0.27)
Linear regression |  |  |  |
on test data | 0.0306 (0.20) | 0.0286 (0.19) | 0.0235 (0.16) | 0.0214 (0.14)
$1$-day offset |  |  |  |
Kalman filter | 0.0279 (0.19) | 0.0278 (0.19) | 0.0247 (0.16) | 0.0211 (0.14)
$3$-day offset |  |  |  |
Kalman filter | 0.0286 (0.19) | 0.0300 (0.20) | 0.0244 (0.16) | 0.0233 (0.16)

Residuals are plotted for estimation of the calibration coefficient for
NRLMSIS-00 $1$ day ahead in Fig. 4. Also plotted for comparison is the
predicted standard error for this quantity, $\sqrt{S^{\prime}}$, which adopted
a near constant value. The mean value of $\sqrt{S^{\prime}}$ was $0.0337\times
10^{-12}$ kg.m-3, approximating the RMS value of the residual $0.0279\times
10^{-12}$ kg.m-3. Similar results were found for each model.

Figure 4: Mean orbit density residuals, $r_{\rho}$, for 1-day offset Kalman
filter calibration of NRLMSIS-00 output (solid line). $\pm\sqrt{S^{\prime}}$,
the predicted standard error (dashed lines).

Temporal variation of the orbit mean density residuals is plotted in Fig. 5
for NRLMSISE-00 and Fig. 6 for the combined output. These plots show RMS
residuals computed over a moving window of $100$ orbits, omitting missing
values. In each case, there are periods where the magnitude of residuals for
the test data linear regression and Kalman filter calibrations sharply
increased. Comparing the plots with Fig. 1, these periods often coincided with
increases in density and relatively high geomagnetic activity (large peaks in
$a_{p}$).

Figure 5: RMS mean orbit density residual values, taken over a moving window
of 100 orbits ($\sim 1$ week). Plotted for NRLMSISE-00 training data linear
regression (blue), test data linear regression (orange), 1-day offset Kalman
filter calibration (yellow), and 3-day offset Kalman filter calibration
estimates (purple). Figure 6: RMS mean orbit density residual values, taken
over a moving window of 100 orbits ($\sim 1$ week). Plotted for combined
training data linear regression (blue), test data linear regression (orange),
1-day offset Kalman filter calibration (yellow), and 3-day offset Kalman
filter calibration estimates (purple).

###  4.5 Principal component data

The calibration technique described in Section 3 was also applied to scores
for the first $6$ density principal components. Fixed calibrations of these
scores from linear regressions were also computed for comparison. Predicted
scores were combined to give density estimates as a function of $u$ over each
orbit for each of the empirical models (driven as before). The resulting RMS
residuals are tabulated in Tab. 2. Temporal variation in density residuals is
plotted in Fig. 7 for predictions based on NRLMSISE-00 and Fig. 8 for those
based on a combination of model predictions. Similar sharp increases in
residuals occurred during the same time periods as those for orbit mean
density predictions.

Table 2: RMS residuals for different density principal component calibration methods. All quantities are shown in $10^{-12}$ kg.m-3 (without brackets) and ratio with the mean value (within brackets). Calibration method | NRLMSISE-00 | JB2008 | DTM-2020 | Combined
---|---|---|---|---
Linear regression |  |  |  |
on training data | 0.0581 (0.39) | 0.0508 (0.34) | 0.0582 (0.39) | 0.0521 (0.35)
Linear regression |  |  |  |
on test data | 0.0448 (0.30) | 0.0434 (0.29) | 0.0390 (0.26) | 0.0376 (0.25)
$1$-day offset |  |  |  |
Kalman filter | 0.0457 (0.31) | 0.0431 (0.29) | 0.0382 (0.26) | 0.0364 (0.24)
$3$-day offset |  |  |  |
Kalman filter | 0.0476 (0.32) | 0.0455 (0.30) | 0.0390 (0.26) | 0.0384 (0.26)
Figure 7: RMS density residual values, taken over a moving window of 100
orbits ($\sim 1$ week). Plotted for NRLMSISE-00 test data linear regression
(blue), training data linear regression (orange), 1-day offset Kalman filter
calibration (yellow), and 3-day offset Kalman filter calibration (purple)
principal component estimates. Figure 8: RMS density residual values, taken
over a moving window of 100 orbits ($\sim 1$ week). Plotted for combined test
data linear regression (blue), training data linear regression (orange), 1-day
offset Kalman filter calibration (yellow), and 3-day offset Kalman filter
calibration (purple) principal component estimates.

##  5 Discussion

The results presented in Tab. 1 and Tab. 2 show that the Kalman filter
calibration method provided substantially more accurate predictions than the
linear regression on training data. Allowing the relationship between between
the thermosphere model output and satellite density measurements to evolve
over time helps account for trends not well represented by the climatological
and aerodynamic models used. This was true for each of the models examined and
their combination estimating the calibration coefficients $1$ and $3$ days
ahead.

The proportional improvement using the Kalman filter calibration method was
greater for predicting orbit mean densities than densities at specific
locations (based on principal component scores). Comparison of Fig. 5 and Fig.
7 and of Fig. 6 and Fig. 8 shows very similar trends in RMS residuals over
time. However, there is an approximate floor which is substantially lower for
the orbit mean density residuals. The greater improvement in orbit mean
density predictions reflects the reduced effect of uncertainty in measurements
and effect of unmodelled processes through averaging. Future work will examine
the relative importance of these factors.

The Kalman filter calibration method had similar accuracy to applying a
constant linear regression across all the test data. Moreover, changes in the
accuracy of their predictions over time were highly correlated. Times when
residuals were large may indicate when thermosphere variations not well-
represented by the empirical models were important.

Changes in the calibration coefficients could have resulted from inaccuracies
in thermosphere or satellite aerodynamic modelling. Where the empirical model
response to geomagnetic activity or changing solar flux was inaccurate
calibration changes would have occurred. Empirical thermosphere models all
have greatly reduced performance during geomagnetic storms[25]. The
combination of any inaccuracy in model local time dependence and nodal
precession would have would have produced longer-term changes. Moreover,
unmodelled changes in satellite gas-surface interaction properties and
thermosphere composition would have caused calibration changes by introducing
systematic errors in the POD density measurements. Investigation of the
magnitude of such effects would be worthwhile but is beyond the scope of the
present study.

The accuracy of this technique may be dependent on where in the solar cycle it
is applied. Localised density predictions have previously been reported to be
less accurate at solar minimum, suggesting that stochastic effects originating
from the lower atmosphere are more prominent at this time [4]. Moreover, the
accuracy of Swarm-C POD thermosphere density measurements is reduced at solar
minimum due to relatively low aerodynamic accelerations compared with other
non-gravitational forces [21]. Other satellite-based thermosphere density
estimation techniques would have similar issues. Therefore, it is likely that
better performance could be obtained applying this technique near solar
maximum conditions.

By considering additional data, a better model of the temporal evolution of
the calibration coefficients may be obtained. This could involve introducing
dependences of $R$ and $\mathbf{M}$ on solar flux and geomagnetic indices. The
former would account for change in the magnitude of unmodelled short timescale
thermosphere variability and the latter for variation in model calibration.
Alternatively, the trivial $\mathbf{F}$ used here could be replaced with some
function of these indices.

The methods presented here could form part of a dual Kalman filter approach to
orbit determination and orbit localised density estimation. In such an
approach one Kalman filter would estimate the orbit state and another would
simultaneously estimate the model density calibration. Kalman filtering itself
is computationally inexpensive, which would allow such a technique to be
implemented aboard a satellite. The approach could be implemented in
conjunction with surrogate models representing mean density or principal
component scores in terms of geomagnetic and solar flux indices, reducing
processing and storage costs associated with density models. The development
of operational Kalman filter density prediction tools for satellites would
improve the fidelity of orbit determination and prediction. In so doing it
would aid accurate and efficient collision avoidance and propulsionless
manoeuvring.

##  6 Conclusions

A Kalman filter method has been developed for calibrating thermosphere density
predictions from empirical models for a particular satellite. This method
substantially improved upon density estimates made using a calibration based
on linear regression over previous data. Combination of estimates from
different models using a best linear unbiased technique further improved
performance. However, the Kalman filter calibration technique suffered from
periods of reduced performance, suggestive of short time-scale thermosphere
processes which were not well-represented by the climatological empirical
models.

## Acknowledgements

This work was carried out with funding from the Royal Australian Air Force.
Resources from the National Computational Infrastructure (NCI), supported by
the Australian Government, were used. Thanks go to Dr. Brenton Smith and Dr.
Li Qiao for helpful discussions regarding astrodynamics and Kalman filtering.
Thanks also to Prof. Russell Boyce and Dr. Melrose Brown for advice and
general guidance on preparation of this paper.

## References

  * [1] D. A. Vallado, D. Finkleman, A critical assessment of satellite drag and atmospheric density modeling, Acta Astronautica 95 (2014) 141–165.
  * [2] J. Emmert, Thermospheric mass density: A review, Advances in Space Research 56 (5) (2015) 773–824.
  * [3] N. B. Stastny, F. R. Chavez, C. Lin, T. A. Lovell, R. A. Bettinger, J. Luck, Localized density/drag prediction for improved onboard orbit propagation, Tech. rep., AIR FORCE RESEARCH LAB KIRTLAND AFB NM SPACE VEHICLES DIRECTORATE (2009).
  * [4] D. Pérez, B. Wohlberg, T. A. Lovell, M. Shoemaker, R. Bevilacqua, Orbit-centered atmospheric density prediction using artificial neural networks, Acta Astronautica 98 (2014) 9–23.
  * [5] D. Pérez, R. Bevilacqua, Neural network based calibration of atmospheric density models, Acta Astronautica 110 (2015) 58–76.
  * [6] E. Friis-Christensen, H. Lühr, D. Knudsen, R. Haagmans, Swarm–an earth observation mission investigating geospace, Advances in Space Research 41 (1) (2008) 210–216.
  * [7] M. Mahooti, NRLMSISE-00 Atmosphere Model, [Online; accessed 8-September-2020] (2019).
URL
https://www.mathworks.com/matlabcentral/fileexchange/56253-nrlmsise-00-atmosphere-
model

  * [8] M. Mahooti, Jacchia-Bowman Atmospheric Density Model, [Online; accessed 30-April-2021] (2018).
URL https://www.mathworks.com/matlabcentral/fileexchange/56163-jacchia-bowman-
atmospheric-density-model

  * [9] D. Dalaur, S. Bruinsma, D. Jackson, SWAMI MCM Model, [Online; accessed 2-July-2021] (2021).
URL https://github.com/swami-h2020-eu/mcm

  * [10] J. Picone, A. Hedin, D. P. Drob, A. Aikin, NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues, Journal of Geophysical Research: Space Physics 107 (A12) (2002) SIA–15.
  * [11] A. Hedin, J. Salah, J. Evans, C. Reber, G. Newton, N. Spencer, D. C. Kayser, D. Alcayde, P. Bauer, L. Cogger, et al., A global thermospheric model based on mass spectrometer and incoherent scatter data MSIS, 1. N2 density and temperature, Journal of Geophysical Research 82 (16) (1977) 2139–2147.
  * [12] A. Hedin, C. Reber, G. Newton, N. Spencer, H. Brinton, H. Mayr, W. Potter, A global thermospheric model based on mass spectrometer and incoherent scatter data MSIS, 2. Composition, Journal of Geophysical Research 82 (16) (1977) 2148–2156.
  * [13] B. Bowman, W. K. Tobiska, F. Marcos, C. Huang, C. Lin, W. Burke, A new empirical thermospheric density model JB2008 using new solar and geomagnetic indices, in: AIAA/AAS astrodynamics specialist conference and exhibit, 2008, p. 6438.
  * [14] S. Bruinsma, The DTM2020 models, in: 43rd COSPAR Scientific Assembly. Held 28 January-4 February, Vol. 43, 2021, p. 813.
  * [15] F. Barlier, C. Berger, J. Falin, G. Kockarts, G. Thuillier, A thermospheric model based on satellite drag data., in: Annales de Geophysique, Vol. 34, 1978, pp. 9–24.
  * [16] S. Bruinsma, The DTM-2013 thermosphere model, Journal of Space Weather and Space Climate 5 (2015) A1.
  * [17] R. E. Kalman, A new approach to linear filtering and prediction problems, Transaction of the ASME-Journal of Basic Engineering (1960) 35–45.
  * [18] R. Mehra, Approaches to adaptive filtering, IEEE Transactions on automatic control 17 (5) (1972) 693–698.
  * [19] J. C. Pinheiro, D. M. Bates, Unconstrained parametrizations for variance-covariance matrices, Statistics and computing 6 (3) (1996) 289–296.
  * [20] L. Lyons, D. Gibaut, P. Clifford, How to combine correlated estimates of a single physical quantity, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 270 (1) (1988) 110–117.
  * [21] J. van den IJssel, E. Doornbos, E. Iorfida, G. March, C. Siemes, O. Montenbruck, Thermosphere densities derived from swarm gps observations, Advances in Space Research 65 (7) (2020) 1758–1771.
  * [22] Collecte Localisation Satellites, Space weather services at CLS (2015).
URL https://spaceweather.cls.fr/

  * [23] GFZ German Research Centre for Geosciences, Geomagnetic Kp Index (2020).
URL https://www.gfz-potsdam.de/en/kp-index/

  * [24] Space Environment Technologies, JB2008 (2021).
URL https://sol.spacenvironment.net/JB2008/index.php

  * [25] S. Bruinsma, C. Boniface, E. K. Sutton, M. Fedrizzi, Thermosphere modeling capabilities assessment: geomagnetic storms, Journal of Space Weather and Space Climate 11 (2021) 12.