# VAE-IF: Deep feature extraction with averaging for unsupervised artifact
detection in routine acquired ICU time-series ††thanks: This research was
funded in part by an EPSRC Doctoral Training Partnership PhD studentship to HH
and The Carnegie Trust for the Universities of Scotland (RIG009251) to JE. For
the purpose of open access, the author has applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

Hollan Haule, Ian Piper, Patricia Jones, Chen Qin, Tsz-Yan Milly Lo, and
Javier Escudero H. Haule and J. Escudero are with the Institute for Imaging,
Data and Communications (IDCOM), School of Engineering, University of
Edinburgh, Edinburgh EH9 3FB, U.K. (e-mail: hhaule@ed.ac.uk). T.-Y.M. Lo and
I. Piper are with the Centre of Medical Informatics, Usher Institute,
University of Edinburgh, U.K., and also with the Royal Hospital for Children
and Young People, Edinburgh, U.K. C. Qin is with Department of Electrical and
Electronics Engineering & I-X, Imperial College London, U.K. P. Jones is with
Department of Child Life and Health, University of Edinburgh, Edinburgh,
U.K.Manuscript submitted December 2023.

###### Abstract

Artifacts are a common problem in physiological time-series data collected
from intensive care units (ICU) and other settings. They affect the quality
and reliability of clinical research and patient care. Manual annotation of
artifacts is costly and time-consuming, rendering it impractical. Automated
methods are desired. Here, we propose a novel unsupervised approach to detect
artifacts in clinical-standard minute-by-minute resolution ICU data without
any prior labeling or signal-specific knowledge. Our approach combines a
variational autoencoder (VAE) and an isolation forest (iForest) model to learn
features and identify anomalies in different types of vital signs, such as
blood pressure, heart rate, and intracranial pressure. We evaluate our
approach on a real-world ICU dataset and compare it with supervised models
based on long short-term memory (LSTM) and XGBoost. We show that our approach
achieves comparable sensitivity and generalizes well to an external dataset.
We also visualize the latent space learned by the VAE and demonstrate its
ability to disentangle clean and noisy samples. Our approach offers a
promising solution for cleaning ICU data in clinical research and practice
without the need for any labels whatsoever.

###### Index Terms:

Unsupervised learning, Autoencoders, Artifact detection, Intensive care unit
data, Isolation forests.

##  I Introduction

Continuous monitoring of vital signs in Intensive Care Units (ICU) is
essential for patient management and treatment. These vitals typically include
signals such as blood pressure, heart rate, and temperature. The clinical
practice in most ICU involves the routine acquisition of these signals at low
temporal resolution (minute-by-minute samples). This data can become a very
useful resource for medical research, as it could advance knowledge about
disease processes and patient care [1]. However, extracting useful insights
from this data is a big challenge due to artifacts.

Routine collection of ICU data is done in uncontrolled conditions. The focus
is clinical care rather than medical research [2], which often affects data
quality. By way of illustration, probe changes and accidental probe
displacement, position changes, and clinical interventions [3] can manifest as
“spikes” or “flat-line” artifacts in the collected data, among other forms of
artifacts. Therefore, re-purposing this data for offline clinical research
risks getting misleading results if artifacts are not dealt with.

The gold standard for dealing with data artifacts is manual annotation by
experienced researchers, but this is costly and non-scalable [4]. Therefore,
developing techniques to help automate this task is an active and important
area of research. There are several approaches for detecting artifacts in
diverse types of physiological time series. Some of these approaches include
Signal Quality Indices (SQI) [5, 6, 7, 8, 9] based on template matching,
frequency-based methods [10], and machine learning (ML) [11, 12, 13, 14, 15,
16, 17]. Unsupervised ML techniques are of special interest in this domain
because of the scarcity of labeled data. However, existing approaches either
are signal-specific or work only on high-temporal resolution, quasi-periodic
physiological signals, which leaves a gap for methods that work on non-
periodic, routinely collected physiologic data from ICU clinical practice.

We aim to develop unsupervised ML methods for automated artifact detection in
minute-by-minute ICU time series. Building upon a labeled real-world dataset,
we initially explore supervised ML techniques known for their reliable
performance with labeled data. Then, we introduce a novel completely
unsupervised ML approach that exhibits comparable performance to the fully
supervised one. We validate the generalizability of this approach on
recordings from an external dataset. Our methods leverage concepts from LSTMs,
variational autoencoders (VAE), and decision trees. In particular, our
contributions are:

  1. 1.

A completely unsupervised approach to detect artifacts in minute-by-minute ICU
data samples based on VAE and isolation Forest model (VAE-IF) that does not
require previous labelling of artifacts and exploits averaged latent sample
representations for enhanced accuracy.

  2. 2.

A supervised model based on LSTM with self-attention (LSTM-ATTn) as an
alternative for artifact detection when labels are available in a training
set.

  3. 3.

Experiments with a real-world ICU dataset to compare the performance of the
alternative models for artifact detection and visualise the latent
representations by the unsupervised approach.

The manuscript is organised as follows. Section II covers literature review.
Section III introduces our methodology and benchmark techniques. Section IV
then presents the experiments and results. Finally, Sections V to VII discuss
and conclude our study.

##  II Literature review

Some of the earliest works on artifact detection and removal in ICU
physiological time series data used a combination of statistical rule-based
methods [18] and wavelet transform [10]. Jakob et al. [18] compared Rosner
statistic, median detection, and a rule-based slope detection method to detect
artifacts in physiological data from a paediatric ICU. On average, the three
methods achieved sensitivity of 66%, 72%, and 64% respectively. The authors
concluded that median filtering is efficient against spikes/outliers but not
suitable for artifacts spanning a long duration. Saeed et al. [10] used the
Haar wavelet transform to detect artifacts and physiological trend changes in
multivariate physiological signals. The method used normalised wavelet
coefficients to detect significant trend changes in the signals. Using scale
4, which represented changes in 160-second intervals, and XOR logic the method
could detect physiological implausible changes as artifacts. 12 out of 15
segments containing artifacts were detected with no false positives. However,
the method relies on trend changes and may make it ineffective against
flatline artifacts.

SQI is among the most prominent techniques encountered in the literature about
physiological time series artifact detection and removal [5, 6, 7, 8, 9]. SQI
estimates signal quality based on the morphology of a pulsatile physiological
signal. SQI helps identify and exclude segments corrupted by noise when
deriving vital signs. Zong et al. [19] developed an SQI combining waveform
morphological features and fuzzy logic to estimate the quality of arterial
blood pressure (ABP) waveform, a continuous value ranging between 0 and 1,
which was then used to estimate and reduce blood pressure false alarms. Sun et
al. [8] developed a binary-valued SQI for ABP combining morphological features
and physiologic constraints, where 0 and 1 indicate clean and noisy beats,
respectively. Apart from SQI for blood pressure signal, Yaghouby et al. [5]
developed an SQI for electrocardiogram signals that uses a correlation
coefficient between a beat detected in a segment and template beat as an
indicator for signal quality. Silva et al. [6] developed an SQI for multi-
channel physiologic recordings using a multi-channel adaptive filter and
applied it to photoplethysmogram (PPG) signal. Though this SQI offers some
flexibility, it is important to note that, in general, SQI are tuned to work
on specific signals and are designed to be highly specialized. More
importantly, they require high-temporal resolution data.

ML techniques have also been used to detect and remove artifacts ICU
physiology data. One of the earliest ML-based work in artifact detection for
ICU minute resolution time series data is Tsien et al. [11]. This study
formulated the task as a supervised learning problem using a decision tree
model and a dataset with annotations by an expert. The time series were
segmented into overlapping segments for which a feature vector was computed.
Labels were assigned to each segment depending on whether more or less than
50% of samples were artifacts. However, this approach potentially wastes data
especially for long window sizes.

More recently, numerous ML techniques for artifact detection in ICU
physiological time-series have been explored [12, 13, 14, 15, 16, 17].
Mataczynski et al. [12] used a deep learning model based on ResNets to
classify Intracranial Pressure (ICP) waveform pulses as: normal, possibly
pathological, likely pathological, pathological or artifact. Subramanian et
al. [13] compared unsupervised ML models: $k$-nearest neighbours, isolation
forests, and 1-class support vector machines (OCSVM) to detect artifacts in
electrodermal activity data sampled at 256Hz. The time series were transformed
into 0.5s windows and 12 features were extracted. The unsupervised models
computed scores reflecting the presence of artifacts in the window. Chen et
al. [14] use a convolutional neural network (CNN) on 10s short-time Fourier
transformed segments of a PPG signal to classify them as good or bad.
Edinburgh et al. [17] used an unsupervised model based on variational
autoencoders (VAE) to detect artifacts in 10s windows of invasive arterial
blood pressure. All of these techniques work on high-resolution data and
detect artifacts at window/segment level. Hence, they are not directly
applicable to routinely collected low-resolution clinical signals.

In contrast, Wu et al. [16] worked on relatively low-resolution invasive blood
pressure time series with 1 sample every 15s and predicted the presence of
artifacts at the sample level. The authors compared forecasting methods
(ARIMA, exponential smoothing), XGBoost, isolation forest, and OCSVM, and
found that XGBoost significantly outperformed other methods with over 99%
sensitivity and specificity. However, XGBoost being a fully-supervised model
requires labeled data to train. In addition, the study showed that
multivariate models which take in mean, systolic and diastolic blood pressure
outperformed univariate models which take in systolic blood pressure only.

It is worth considering techniques developed for detecting anomalies in other
kinds of time series as well since this problem closely resembles artifact
detection in physiological time series. Munir et al. [20] developed a CNN-
based model with a predictor and anomaly detector module. The predictor module
predicts the next timestamp sample, which is then fed into the anomaly
detector module. The anomaly detector module computed the Euclidian distance
between actual sample and predicted sample, which determined if a sample is an
outlier based on a set threshold. The method could detect both point and
contextual anomalies in periodic time-series.

In this work, we develop a completely unsupervised approach utilizing
$\beta$-VAE [21] and Isolation forest [22] to detect artifacts at sample level
in minute-resolution vital signs data, including mean blood pressure (BPm),
heart rate (HRT) and mean intracranial pressure (ICPm), aiming to completely
eliminate the need for labeled data when building models for artifact
detection. Unlike previous studies utilizing VAEs to detect
artifacts/anomalies by first training on clean data, our method trains on raw
data containing both clean samples and artifacts. We also develop a fully
supervised classifier based on LSTM with attention as a proposed benchmark
[23, 24]. We also implement isolation forest on its own and XGBoost as
benchmarks. We perform an ablation study to demonstrate the significance of
feature extraction with a VAE and the importance of the self-attention
component in the LSTM classifier. We also train an LSTM classier without the
self-attention component to demonstrate its significance. Additionally, we
explore the latent space of real-world data containing artifacts. In contrast
with previous techniques, we train our model on a dataset of all the mentioned
vital signs mixed. We show that the model learns and generalizes the dynamics
of artifacts across the different vital signs. We validate our method using
gold-standard annotations by an experienced expert and demonstrate the
generalizability of the method on the MIMIC-IV [25] dataset without
retraining.

##  III Methodology

In this Section, we present a novel unsupervised approach for artifact
detection in minute-by-minute physiological ICU recordings, addressing the
challenge of limited labeled data. As highlighted in Section I, the scarcity
of labeled data makes our approach particularly valuable. While we recognize
the scarcity of labeled data, we first establish a performance benchmark using
supervised models. Then, we present and evaluate our novel unsupervised
approach, showcasing its generalizability to an external dataset. However,
before delving into the details of the models, we describe the dataset,
preprocessing, and problem formulation.

###  III-A KidsBrainIT dataset

We have access to a research dataset from the KidsBrainIT project [26] that
was collected prospectively and extracted for research purposes from Paeditric
Intensive Care Units (PICU). We use the KidsBrainIT pilot dataset, which
contains about 98 patients from 2 PICU in the UK and about 27 patients that
were added later on, making up a total of 125. Each patient in the dataset has
around 14 physiological signals that were collected at minute-by-minute
intervals, a routine clinical setup in ICU. It also includes annotations on
artifacts by an expert researcher. For our experiments, we used the BPm, ICPm,
and HRT as these signals have the least amount of missing values across all
patients recordings. Fig. 1 shows the distribution of artifacts and good
signal for the selected vital signs.

Figure 1: Distribution of amplitude values for artifactual and non-artifactual
samples in KidsBrainIT dataset.

###  III-B Data preprocessing

We selected 64 patients out of the 125, with at least 90% data in BPm, ICPm,
and HRT recordings. We then randomly split these patients into 47 and 17 as
training and testing sets respectively. We standardised all patient recordings
individually by subtracting the median followed by division with the
Interquartile range of the recording. Standardizing each recording
individually ensures that the differences in vital sign dynamics between
patients are preserved. We also applied quadratic interpolation to fill in the
missing values in the standardized recordings. We created a sequence of sample
labels from each patient recording based on annotations by the expert, i.e. 1
for “artifact” and 0 otherwise. We use these labels as ground truth to measure
the performance of different models. To preserve privacy, we shifted the
actual timestamps of the recordings to midnight 1 Jan, 2010. Fig. 2 shows an
example of a pre-processed patient recording. We then applied an overlapping
sliding window of size $W$ minutes, and step size of 1 minute (1 sample) on
the recordings across all patients in the training set. $W$ is a
hyperparameter, that is tuned during the experiments.

Figure 2: Sample patient standardised recording, containing three signals
(BPm, ICPm, HRT) and their corresponding artifact labels, (BPm_label,
ICP_label, HRT_label), as annotated by an expert researcher.

###  III-C Problem formulation

We consider a preprocessed dataset of size $N$,
$D:=\left\\{\left(\textbf{x}_{i},\textbf{y}_{i}\right)\right\\}$ such that
$\textbf{x}\in\mathbb{R}^{W},\textbf{y}\in\left\\{0,1\right\\}^{W},\mathit{i}\in\left\\{1\ldots
N\right\\}$, where $\textbf{x}_{i}$ is a segment of length $W$ of the time-
series and $\textbf{y}_{i}$ is the corresponding label for each sample in the
segment as described in Section III-B. Additionally, we define
$\textbf{z}\in\mathbb{R}^{H}$, a vector of size $H$, as the latent
representation of a sample $x$ in the segment x. The latent representation for
a signal of length $L$ is represented by a matrix
$\textbf{Z}\in\mathbb{R}^{L\times H}$.

###  III-D Proposed benchmarks

In this section, we describe the supervised models to be used as a performance
benchmarks for the unsupervised approach. For this task, we build two models:
a neural-network model based on LSTM with self-attention, and XGBoost.

####  III-D1 LSTM with self-attention (LSTM-ATTn)

LSTM [23] is a type of recurrent neural network architecture for learning
sequential data. An LSTM can learn temporal information over long sequences
through a hidden state. The hidden state at a time step is given by the
following equation,

$\mathbf{h}_{t}=LSTM(x_{t},\mathbf{h}_{t-1}),$ (1)

where $\mathbf{h}_{t}\in\mathbb{R}^{H}$ is the hidden state vector and
$x_{i}\in\mathbb{R}$ is a sample in segment $\mathbf{x}$ as defined in III-C.
A Bidirectional LSTM [27] has access to past and future information. Its
hidden state has an effective size of $\mathbf{h}_{t}\in\mathbb{R}^{2H}$.

We stack two LSTM layers with hidden size $H$ followed by a scaled dot-product
attention [24], given by

${\mathbf{h}}^{\prime}_{i}=\sum_{j}^{W}\mathrm{softmax}\left(\frac{\mathbf{h}_{i}^{T}\mathbf{h}_{j}}{\sqrt{H}}\right)\mathbf{h}_{j},$
(2)

where ${\mathbf{h}}^{\prime}_{i}$ is the final hidden state and $W$ is the
window size as defined in Section III-B. We also explore the effect of using a
bidirectional LSTM during hyper-parameter tuning. Next, we add a dense layer
with a sigmoid activation function to give a probability estimate
${y}^{\prime}_{i}$ of sample $x_{i}$ being an artifact. We train the model
using Adam optimizer [28]. Fig 3 shows the architecture diagram of the LSTM-
Attn classifier.

####  III-D2 XGBoost

XGBoost [29], or Extreme Gradient Boosting, is a decision tree ensemble based
on CART (Classification and Regression Trees) models trained using the
gradient boosting algorithm [30]. The algorithm trains additively by fitting
new trees to the residuals. XGBoost is known for its excellent performance
[16], and we use it as a second benchmark.

###  III-E Fully unsupervised approach: VAE-IF

In this section, we present our proposed unsupervised approach for artifact
detection in minute-by-minute physiological recordings. Notably, this approach
is particularly valuable as it does not rely on labeled data. At a high level,
our approach involves pre-training a VAE [31] model on the preprocessed
training dataset, which includes both noisy and clean segments,
$\mathbf{x}_{i}$, as described in Sections III-B and III-C. Subsequently, we
utilize the encoder of the trained VAE to extract features from a test signal.
Once the features are extracted, we directly employ the Isolation Forest [22]
anomaly detection algorithm to classify samples in the signal as normal or
artifact (anomaly). In the following subsections, we give a detailed
description of the main components of this approach.

####  III-E1 Training the VAE

Our model is based on $\beta$-VAE [21], a variant of VAE with the following
loss function (Evidence Lower Bound, ELBO) in Eq. 3:

$\mathcal{L}(\theta,\phi;\mathbf{x})=\mathbb{E}_{q_{\phi}(\mathbf{z|x})}\left[\log
p_{\theta}(\mathbf{x|z})\right]-\beta
D_{KL}(q_{\phi}(\mathbf{z|x})||p_{\theta}(\mathbf{z})),$ (3)

where $\theta,\phi$ are the parameters of the decoder
$p(\mathbf{x}|\mathbf{z})$ and encoder $q(\mathbf{z}|\mathbf{x})$ and $\beta$
is a parameter that can be tuned to encourage the model to learn disentangled
representations in the latent space. The encoder consists of two stacked LSTM
layers, followed by attention and a dense layer. Similarly, the decoder has a
comparable architecture, except for the absence of the attention block. The
first three layers of the encoder are similar to that of the LSTM-Attn model
in Fig 3. The full architecture for the VAE can be seen in Fig. 4.

We train the VAE with a mix of clean and artifact-containing data segments x.
We use Bayesian optimization (BO) [32] to tune model hyper-parameters listed
in Table I.

Figure 3: Architectual diagram of the supervised LSTM-Attn benchmark Figure 4:
Architectural diagram of the VAE

####  III-E2 Feature extraction and averaging

We use the encoder of a trained VAE to extract features from a target
physiological time series. We apply a sliding window with unit step size on
the target signal, resulting in overlapping segments. Passing these
overlapping segments through the encoder to extract features will result in
some samples having more than one embedding vector. This is crucial since
averaging these embeddings gives a more accurate representation of the sample.
This operation is summarized in Fig. 5. Essentially, each sample in the
original signal is converted into a vector that represents the state of the
sample given the previous samples. The transformed sample windows are then
averaged to obtain a sequence of state vectors $\bar{Z}$.

$x_{1}$$x_{2}$…$x_{w}$$\acute{x}_{2}$$x_{3}$…${x}_{w+1}$$z_{1}^{H}$$z_{2}^{H}$…$z_{w}^{H}$$\acute{z}_{2}^{H}$$\acute{z}_{3}^{H}$…$\acute{z}_{w+1}^{H}$$z_{1}^{2}$$z_{2}^{2}$…$z_{w}^{2}$$\acute{z}_{2}^{2}$$\acute{z}_{3}^{2}$…$\acute{z}_{w+1}^{2}$$z_{1}^{1}$$z_{2}^{1}$…$z_{w}^{1}$$\acute{z}_{2}^{1}$$\acute{z}_{3}^{1}$…$\acute{z}_{w+1}^{1}$EncoderFeature
Extraction$\acute{z}_{2}^{H}$…$\acute{z}_{w}^{H}$⋮…⋮$\acute{z}_{2}^{1}$…$\acute{z}_{w}^{1}$$\acute{z}_{2}^{1}$…$\acute{z}_{w}^{1}$+$z_{1}^{H}$$z_{2}^{H}$…$z_{w}^{H}$$z_{w+1}^{H}$⋮⋮…⋮⋮$z_{1}^{2}$$z_{2}^{2}$…$z_{w}^{2}$$z_{w+1}^{2}$$z_{1}^{1}$$z_{2}^{1}$…$z_{w}^{1}$$z_{w+1}^{1}$$\bar{Z}$Isolation
ForestArtifact detectionArtefact01 Figure 5: Feature extraction via
overlapping sliding-window on target signal $x$. This process employs a
trained VAE for encoding, leading to multiple latent representations (e.g.,
$z$ and $z^{\prime}$) for some signal samples based on their position. These
representations are averaged to produce a consistent signal representation,
$\bar{Z}$, which is subsequently input to the iForest model for artifact
detection.

####  III-E3 Artifact detection with iForest

Isolation Forest (iForest) [22] is an algorithm for anomaly detection which
uses an ensemble of isolation trees to partition data based on the assumption
that anomalies exhibit distinct characteristics that set them apart from
normal data points. Consequently, anomalies are expected to be isolated closer
to the root node of the trees compared to normal data points. Compared to
other algorithms like OCSVM, iForest is computationally faster and scales well
to high-dimensional datasets [22].

We use Sklearn’s [33] implementation of iForest with its default
hyperparameters for artifact detection. iForest produce a sequence with output
values corresponding to clean and artifactual samples, represented by the
vector sequence $\bar{Z}$. We first train the iForest model on the extracted
features of all patients and the corresponding target signals in training set
before using it on the testing set.

TABLE I: Hyperparameters and ranges included in the BO

(a) VAE-IF Hyperparameter | Range
---|---
lr | [1e-4, 0.001] (log scale)
latent_dim | [2, 90]
$\beta$ | [1e-3, 1e2] (log scale)
n_epochs | [8, 30]
window_size | [15, 60]
bidirectional | {True, False}

(b) LSTM-ATTn Hyperparameter | Range
---|---
lr | [1e-4, 0.001] (log scale)
n_epochs | [8, 30]
window_size | [15, 60]
bidirectional | {True, False}

(c) XGBoost Hyperparameter | Range
---|---
lr | [0.01, 0.5]
n_estimators | [100, 1000]
window_size | [15, 60]

###  III-F Ablation study

We also perform an ablation study to investigate the individual contributions
and significance of various components in our proposed methods. By
systematically removing specific elements and evaluating their impact on the
performance, we gain insights into the effectiveness of each component.

####  III-F1 Feature extraction ablation

To assess the importance of feature extraction using $\beta$-VAE [21], we
compare the performance of model with a vanilla Isolation Forest model trained
directly on the patient recordings. By excluding the VAE-based feature
extraction step, we examine whether the extracted features significantly
contribute to the detection of artifacts in minute-resolution vital signs
data.

####  III-F2 Self-Attention Ablation

We also analyze the significance of the self-attention component in the
encoder of the VAE. We exclude the self-attention component in the encoder of
the VAE and train the model with the same set of hyper-parameters estimated by
BO. This ablation study helps us understand the added benefit of incorporating
self-attention for accurate artifact detection.

##  IV Experiments and Results

We used BO to tune model hyperparameters for the VAE part of the VAE-IF, LSTM-
ATTn, and XGBoost, using Geometric Mean as the objective function. We also
explored training hyperparameters such as window size ($W$), learning rate,
and number of epochs. The range of hyperparameters explored for the three
models is summarised in Table I.

Using the test set, we evaluate the performance of VAE-IF and compare it to
the supervised models; LSTM-ATTn and XGBoost. To measure performance, we
compute sensitivity and specificity on a sample-by-sample basis. We compute
the performance metrics from each recording in the test set separately and
report the mean and standard deviation scored by each model.

We visualize the latent space learned by the VAE using t-SNE [34] to show how
clean signals and artifactual signals are represented in the latent space. We
use default parameters for t-SNE, with the exception of perplexity which we
set to 5.

We use Pytorch and Pyro [35] to implement the VAE. BO was performed using
ax.dev [36], an open-source platform for managing adaptive experiments. We
trained our models on Nvidia’s GeForce RTX 3080 GPUs.

###  IV-A Unsupervised approach

We performed BO to tune hyperparameters listed in Table I. Tuning
hyperparameters was done in 60 iterations, whereby in each iteration, 70$\%$
of the training set was used for training and 30$\%$ for validation. The best-
performing hyperparameters were then used to train the VAE and the final
evaluation was done on the test set. The best-performing hyperparameters can
be seen in the Supplementary materials. Table II shows the performance results
of VAE-IF on the test set. The table shows that VAE-IF achieves the highest
sensitivity in two of the three signals.

###  IV-B Benchmarks

Similarly, we tune the hyperparameters for the LSTM-ATTn and XGBoost using BO.
The best hyperparameters for the benchmark models also included in the
supplementary material. Table II and II show the best performance on test set
for LSTM-ATTn and XGBoost, respectively. Results show that LSTM-ATTn achieve
the highest sensitivity for HRT signal. The benchmark models achieved higher
specificity compared to VAE-IF across all three signals.

###  IV-C External validation

We tested the generalizability of VAE-IF model on sample patients from the
MIMIC-IV [25] waveform dataset. Here, we present results for the ICP signal of
patient p15649186. Given the difference in sampling rates between the
datasets, we downsampled the ICP signal from 250Hz to 0.0167Hz (one sample per
minute). We then applied the VAE-IF without any retraining. Fig. 6 shows the
model’s artifact prediction and the patient recording. Additional examples are
included in the supplementary material.

Figure 6: VAE_IF prediction on patient p15649186 from MIMIC-IV waveform
dataset, showing automatic detection of spikes and flatline artifacts without
any retraining.

###  IV-D Latent space visualization with t-SNE

We performed an exploration of the latent space embeddings z based on one of
the patients in KidsBrainIT (N004) by projection to 2D with t-SNE. Using
$W=15$, we first extract embeddings as illustrated in Fig. 5. Each embedding
$\textbf{z}_{i}$ is given a numeric label of value between 0 and 15 (the
optimal window length) depending on how many samples in its corresponding
segment $\textbf{x}_{i}$ are artifacts. We visualise the embeddings in 2D
using t-SNE, shown in Fig. 7.

Figure 7: t-SNE visualization of a sample recording, showing disentangled
latent representation of clean data from artifacts. TABLE II: Artifact
Detection Performance on the KidsBrainIT dataset. Best performance for each
signal is highlighted in bold

(a) VAE-IF Signal | Metric | Mean | Std
---|---|---|---
BPm | Sensitivity | 0.802 | 0.219
| Specificity | 0.921 | 0.041
ICPm | Sensitivity | 0.700 | 0.388
| Specificity | 0.926 | 0.052
HRT | Sensitivity | 0.549 | 0.425
| Specificity | 0.855 | 0.227

(b) LSTM-ATTn Signal | Metric | Mean | Std
---|---|---|---
BPm | Sensitivity | 0.767 | 0.175
| Specificity | 0.991 | 0.010
ICPm | Sensitivity | 0.695 | 0.361
| Specificity | 0.993 | 0.007
HRT | Sensitivity | 0.592 | 0.400
| Specificity | 0.916 | 0.243

(c) XGBoost Signal | Metric | Mean | Std
---|---|---|---
BPm | Sensitivity | 0.774 | 0.163
| Specificity | 0.995 | 0.007
ICPm | Sensitivity | 0.672 | 0.355
| Specificity | 0.990 | 0.015
HRT | Sensitivity | 0.465 | 0.414
| Specificity | 0.915 | 0.210

###  IV-E Ablation study

We first evaluated the importance of feature extraction using the VAE. For
this, we compared the performance of our model with a standard Isolation
Forest model trained directly on the patient recordings. Before training, we
transformed each sample into a feature vector containing the current sample
and 14 previous samples, forming a feature vector of size 15 (i.e., $W=15$ as
was used to train the VAE). Next, we analyzed the significance of the self-
attention component in the encoder of the VAE. To do so, we removed the self-
attention component and trained the model with the same set of hyperparameters
estimated by BO. The results of these ablation studies are presented in Table
III.

TABLE III: Performance Comparison of iForest and VAE-IF without attention, showing that both models achieve low performance compared with the VAE-IF with attention approach.  | iForest | VAE-IF (no attention)
---|---|---
Signal | Metric | Mean | Std | Mean | Std
BPm | Sensitivity | 0.294 | 0.221 | 0.602 | 0.227
| Specificity | 0.973 | 0.020 | 0.920 | 0.040
ICPm | Sensitivity | 0.274 | 0.321 | 0.545 | 0.358
| Specificity | 0.959 | 0.030 | 0.937 | 0.042
HRT | Sensitivity | 0.270 | 0.366 | 0.395 | 0.391
| Specificity | 0.885 | 0.255 | 0.863 | 0.206

##  V Discussion

In this study, we proposed an unsupervised approach leveraging the power of
VAE and iIForest for the detection of artifacts in physiological time series
collected at minute-by-minute resolution in ICU settings. To assess its
effectiveness, we conducted a comparative analysis against LSTM-Attention
(LSTM-ATTn) and XGBoost, both of which are supervised models. Our results
demonstrate that our proposed VAE-IF method achieves sensitivity levels
comparable to LSTM-ATTn and XGBoost while maintaining high specificity.
Notably, it can operate in circumstances where no labels about artifacts are
available. Our results clearly demonstrate the ability of VAE-IF to accurately
detect artifacts. Furthermore, we visualize the learned latent space
representation using t-SNE, revealing that a VAE trained on raw data
containing both clean and artifactual samples effectively disentangles clean
from noisy instances. These findings highlight the potential of unsupervised
methods for artifact detection in physiological time series data, offering a
promising solution to the scarcity of labeled data.

In our approach, we employ feature extraction with averaging. Understanding
that the latent state of a sample might be influenced by its neighboring
samples and their relative positions, we utilized an overlapping sliding-
window. When combined with accessing both past and future states via
bidirectional LSTM, averaging leads to a more accurate latent representation
by taking into account multiple neighborhood contexts.

In comparison to the existing literature [11, 12, 13, 14, 17], which detects
artifacts at segment/window level, our VAE-IF approach detects artifacts on a
sample-by-sample basis. It is, however, critical to acknowledge that most of
these studies focused on waveform data. Here, the primary task involves the
identification of deformed morphology in a periodic waveform pulse. In such
contexts, it is logical to filter out the entire segment if it is deformed. By
contrast, our sample-by-sample strategy is particularly relevant for low-
resolution, non-periodic signals. Nonetheless, our approach achieves
sensitivity levels similar to XGBoost, a model known for achieving high
performance as demonstrated in a previous study [16], which also used a
sample-by-sample strategy.

Our approach contributes to the existing literature by presenting a compelling
case training a VAE on raw data comprising both normal (clean) and abnormal
(artifacts) data. In the domain of unsupervised anomaly detection using auto-
encoders, which is closely related to artifact detection, the conventional
practice entails training the model solely on clean data [37, 38] and
subsequently detecting anomalies during testing by measuring the deviation
from the normal data using a distance metric. However, this training strategy
becomes counterproductive when labeled data is scarce. In contrast, our
approach offers a more practical solution for cleaning ICU data in clinical
research. By training the model directly on raw data, we avoid the need for
labeled data while still achieving comparable performance. This highlights the
potential of our approach as an effective tool for cleaning ICU data, enabling
more efficient and reliable clinical research.

In terms of practical utility, multiple factors contribute to the
effectiveness of a model for artifact detection. These factors encompass
performance, ease of model training, and availability of training examples.
Notably, all the models developed in this study exhibit fast training times
when utilizing a GPU. Regarding performance, our unsupervised approach
demonstrates sensitivity levels comparable to other models, albeit with
slightly lower specificity. An additional advantage of our approach is that it
does not necessitate any labeled data for training, enhancing its practical
applicability.

##  VI Limitations and future work

One noteworthy finding of this study is that all models achieved lower
sensitivity values when detecting artifacts in the ICPm and HRT signals
compared to the BPm signal. A possible explanation for this observation is
that a considerable proportion of artifacts in the ICPm and HRT signals share
values within the range of normal data, as illustrated in Fig. 1. Fig. 2
presents an illustrative example of an ICPm recording with an artifact
exhibiting the same value range as normal data, illustrating the challenges
associated with accurately detecting such artifacts.

The results of our study reveal a potential limitation of our approach, as it
demonstrates a slightly lower specificity compared to the supervised models.
However, it is important to consider a counter-argument in the context of
artifact detection. This lower specificity may still be deemed acceptable, as
it could indicate the presence of missed artifacts during the manual
annotation process. It highlights the possibility of capturing additional
instances that human annotators might have overlooked. Another limitation of
our approach is that the model cannot be trained end-to-end, which may hinder
the potential for improved performance. Nevertheless, this limitation presents
an interesting avenue for future research, specifically in the development of
a “deep” isolation forest model tailored specifically for artifact detection.
Exploring such a direction could potentially enhance the performance and
effectiveness of our approach in identifying artifacts with greater accuracy
and reliability.

##  VII Conclusion

We have successfully developed an unsupervised approach for artifact detection
in minute-by-minute physiological recordings, utilizing real-world, routinely
collected ICU data. Our approach has demonstrated performance on par with
fully supervised methods, highlighting its efficacy and potential for
practical implementation as it does not require any data labelling.
Furthermore, we have validated the generalizability of our approach by testing
it on sample patients from an external dataset, and the model predictions
align closely with visual observations. Future studies involving diverse
datasets will provide valuable insights and ensure the robustness of our
approach in various clinical settings. Overall, our findings present a
promising direction for advancing artifact detection methodologies in
physiological time series data, with potential applications in clinical
research.

## Acknowledgment

The authors thank Rob Donald for the discussions on artifact detection using
the KidsBrainIT-pilot dataset.

## References

  * [1] B. Depreitere, F. Güiza, G. Van den Berghe, M. U. Schuhmann, G. Maier, I. Piper, and G. Meyfroidt, “Can optimal cerebral perfusion pressure in patients with severe traumatic brain injury be calculated based on minute-by-minute data monitoring?” _Acta Neurochirurgica, Supplementum_ , vol. 122, pp. 245–248, 2016.
  * [2] D. Shillan, J. A. Sterne, A. Champneys, and B. Gibbison, “Use of machine learning to analyse routinely collected intensive care unit data: A systematic review,” _Critical Care_ , vol. 23, no. 1, 2019.
  * [3] A. E. Johnson, M. M. Ghassemi, S. Nemati, K. E. Niehaus, D. Clifton, and G. D. Clifford, “Machine Learning and Decision Support in Critical Care,” _Proceedings of the IEEE_ , vol. 104, no. 2, 2016.
  * [4] A. Alkhachroum, K. Terilli, M. Megjhani, and S. Park, “Harnessing Big Data in Neurocritical Care in the Era of Precision Medicine,” 2020.
  * [5] F. Yaghouby, C. Daluwatte, and C. G. Scully, “Variability analysis for noisy physiological signals: A simulation study,” in _2016 IEEE Signal Processing in Medicine and Biology Symposium, SPMB 2016 - Proceedings_ , 2017.
  * [6] I. Silva, J. Lee, and R. Mark, “Photoplethysmograph quality estimation through multichannel filtering,” in _Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBS_ , 2011.
  * [7] Q. Li and G. D. Clifford, “Suppress false arrhythmia alarms of ICU monitors using heart rate estimation based on combined arterial blood pressure and ECG analysis,” in _2nd International Conference on Bioinformatics and Biomedical Engineering, iCBBE 2008_ , 2008, pp. 2185–2187.
  * [8] J. X. Sun, A. T. Reisner, and R. G. Mark, “A signal abnormality index for arterial blood pressure waveforms,” in _Computers in Cardiology_ , vol. 33, 2006, pp. 13–16.
  * [9] W. Zong, G. B. Moody, and R. G. Mark, “Reduction of false arterial blood pressure alarms using signal quality assessment and relationships between the electrocardiogram and arterial blood pressure,” _Medical and Biological Engineering and Computing_ , vol. 42, no. 5, 2004.
  * [10] M. Saeed and R. G. Mark, “Multiparameter trend monitoring and intelligent displays using wavelet analysis,” _Computers in Cardiology_ , pp. 797–800, 2000.
  * [11] C. L. Tsien, I. S. Kohane, and N. McIntosh, “Building ICU artifact detection models with more data in less time.” _Proceedings / AMIA … Annual Symposium. AMIA Symposium_ , pp. 706–710, 2001.
  * [12] C. Mataczynski, A. Kazimierska, A. Uryga, M. Burzynska, A. Rusiecki, and M. Kasprowicz, “End-to-End Automatic Morphological Classification of Intracranial Pressure Pulse Waveforms Using Deep Learning,” _IEEE Journal of Biomedical and Health Informatics_ , vol. 26, no. 2, 2022.
  * [13] S. Subramanian, B. Tseng, R. Barbieri, and E. N. Brown, “Unsupervised Machine Learning Methods for Artifact Removal in Electrodermal Activity,” in _Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBS_ , 2021.
  * [14] J. Chen, K. Sun, Y. Sun, and X. Li, “Signal Quality Assessment of PPG Signals using STFT Time-Frequency Spectra and Deep Learning Approaches,” in _Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBS_ , 2021.
  * [15] S. B. Lee, H. Kim, Y. T. Kim, F. A. Zeiler, P. Smielewski, M. Czosnyka, and D. J. Kim, “Artifact removal from neurophysiological signals: Impact on intracranial and arterial pressure monitoring in traumatic brain injury,” _Journal of Neurosurgery_ , vol. 132, no. 6, 2020.
  * [16] M. Wu, P. Branco, J. X. C. Ke, and D. B. MacDonald, “Artifact Detection in Invasive Blood Pressure Data using Forecasting Methods and Machine Learning,” in _Proceedings - 2020 IEEE International Conference on Bioinformatics and Biomedicine, BIBM 2020_ , 2020.
  * [17] T. Edinburgh, P. Smielewski, S. Eglen, M. Czosnyka, and A. Ercole, “DeepClean \- self-supervised artefact rejection for intensive care waveform data using generative deep learning,” 2019.
  * [18] S. Jakob, I. Korhonen, E. Ruokonen, T. Virtanen, A. Kogan, and J. Takala, “Detection of artifacts in monitored trends in intensive care,” _Computer Methods and Programs in Biomedicine_ , vol. 63, no. 3, pp. 203–209, 2000.
  * [19] W. Zong, G. B. Moody, and R. G. Mark, “Reduction of false arterial blood pressure alarms using signal quality assessment and relationships between the electrocardiogram and arterial blood pressure,” _Medical and Biological Engineering and Computing_ , vol. 42, no. 5, pp. 698–706, 2004.
  * [20] M. Munir, S. A. Siddiqui, A. Dengel, and S. Ahmed, “DeepAnT: A Deep Learning Approach for Unsupervised Anomaly Detection in Time Series,” _IEEE Access_ , vol. 7, 2019.
  * [21] I. Higgins, L. Matthey, A. Pal, C. Burgess, X. Glorot, M. Botvinick, S. Mohamed, and A. Lerchner, “B-VAE: Learning basic visual concepts with a constrained variational framework,” in _5th International Conference on Learning Representations, ICLR 2017_ , 2017.
  * [22] F. T. Liu, K. M. Ting, and Z.-H. Zhou, “Isolation Forest,” in _2008 Eighth IEEE International Conference on Data Mining_ , 2008, pp. 413–422.
  * [23] S. Hochreiter and J. Schmidhuber, “Long Short-Term Memory,” _Neural Computation_ , vol. 9, no. 8, 1997.
  * [24] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Å. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Advances in Neural Information Processing Systems_ , vol. 2017-December, 2017\.
  * [25] B. Moody, S. Hao, B. Gow, T. Pollard, W. Zong, and M. Roger, “MIMIC-IV Waveform Database.”
  * [26] T. Lo, I. Piper, B. Depreitere, G. Meyfroidt, M. Poca, J. Sahuquillo, T. Durduran, P. Enblad, P. Nilsson, A. Ragauskas, K. Kiening, K. Morris, R. Agbeko, R. Levin, J. Weitz, C. Park, and P. Davis, “KidsBrainIT: A new multi-centre, multi-disciplinary, multi-national paediatric brain monitoring collaboration,” in _Acta Neurochirurgica, Supplementum_ , 2018, vol. 126\.
  * [27] M. Schuster and K. K. Paliwal, “Bidirectional recurrent neural networks,” _IEEE Transactions on Signal Processing_ , vol. 45, no. 11, 1997.
  * [28] D. P. Kingma and J. L. Ba, “Adam: A method for stochastic optimization,” in _3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings_ , 2015.
  * [29] T. Chen and C. Guestrin, “XGBoost: A scalable tree boosting system,” in _Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , vol. 13-17-August-2016, 2016.
  * [30] J. H. Friedman, “Greedy function approximation: A gradient boosting machine,” _Annals of Statistics_ , vol. 29, no. 5, 2001.
  * [31] D. P. Kingma and M. Welling, “Auto-Encoding Variational Bayes,” _2nd International Conference on Learning Representations, ICLR 2014_ , 12 2013.
  * [32] J. Mockus, V. Tiesis, and A. Zilinskas, “The application of Bayesian methods for seeking the extremum,” _Towards Global Optimisation_ , vol. 2, no. 2, 1978.
  * [33] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and Ã. Duchesnay, “Scikit-learn: Machine learning in Python,” _Journal of Machine Learning Research_ , vol. 12, 2011.
  * [34] L. Van Der Maaten and G. Hinton, “Visualizing data using t-SNE,” _Journal of Machine Learning Research_ , vol. 9, 2008.
  * [35] E. Bingham, J. P. Chen, M. Jankowiak, F. Obermeyer, N. Pradhan, T. Karaletsos, R. Singh, P. Szerlip, P. Horsfall, and N. D. Goodman, “Pyro: Deep universal probabilistic programming,” _Journal of Machine Learning Research_ , vol. 20, 2019.
  * [36] E. Bakshy, L. Dworkin, B. Karrer, K. Kashin, B. Letham, A. Murthy, and S. S. Facebook, “AE: A domain-agnostic platform for adaptive experimentation,” Tech. Rep., 2018.
  * [37] Z. Xiong, D. Zhu, D. Liu, S. He, and L. Zhao, “Anomaly Detection of Metallurgical Energy Data Based on iForest-AE,” _Applied Sciences (Switzerland)_ , vol. 12, no. 19, 2022.
  * [38] N. Merrill and A. Eskandarian, “Modified Autoencoder Training and Scoring for Robust Unsupervised Anomaly Detection in Deep Learning,” _IEEE Access_ , vol. 8, 2020.


# Fully Three-dimensional Simulation and Modeling of a Dense Plasma Focus

B. T. Meehan, J. H. J. Niederhaus B. T. Meehan is with National Security
Technologies, LLC, a Department of Energy Contractor, e-mail:
<EMAIL_ADDRESS>H. J. Niederhaus is a Computer Scientist at Sandia
National Laboratories, email<EMAIL_ADDRESS>received December
1, 2012; revised January 11, 2013.

###### Abstract

A Dense Plasma Focus (DPF) is a pulsed-power machine that electromagnetically
accelerates and cylindrically compresses a shocked plasma in a Z-pinch. The
pinch results in a brief ($\sim$100 nanosecond) pulse of X-rays, and, for some
working gases, also a pulse of neutrons. A great deal of experimental research
has been done into the physics of DPF reactions, and there exist mathematical
models describing its behavior during the different time phases of the
reaction. Two of the phases, known as the inverse pinch and the rundown, are
approximately governed by magnetohydrodynamics, and there are a number of
well-established codes for simulating these phases in two dimensions or in
three dimensions under the assumption of axial symmetry. There has been little
success, however, in developing fully three-dimensional simulations. In this
work we present three-dimensional simulations of DPF reactions and demonstrate
that 3D simulations predict qualitatively and quantitatively different
behavior than their 2D counterparts. One of the most important quantities to
predict is the time duration between the formation of the gas shock and
Z-pinch, and the 3D simulations more faithfully represent experimental results
for this time duration and are essential for accurate prediction of future
experiments.

###### Index Terms:

Dense Plasma Focus, Magnetohydrodynamics, Simulation and Modeling, Controlled
Fusion

##  I Introduction

There is an extensive literature on experiments performed with Dense Plasma
Focus (DPF) machines, exploring both fundamental Z-pinch physics [1, 2, 3] and
applications of the DPF to fields like X-ray lithography [4], fusion energy
[5], and modeling of astrophysics [6], to name a few. The vast majority of
research has been driven by theory and experimentation, but there is a shift
towards developing new experiments informed by computational models and
simulations. Many of the simulations that have been used to design DPF
experiments are 1-dimensional, in the sense that current, temperature, or
expected neutron yield are computed as a function of a single parameter being
varied [7]. There has also been work to develop codes that are capable of full
magnetohydrodynamic (MHD) modeling of the inverse pinch and run-down phases of
the DPF reaction on spatial domains (see Section II for descriptions of the
reaction phases), and there exist 2D simulations modeling the physics of the
DPF, which are often extended to 3D assuming axial symmetry [8]. Despite some
success with 1D and 2D simulations, fully 3D simulations of the MHD phases of
the DPF – which do not impose symmetry on the physics of the reaction – are
difficult, and the literature presenting the results of such simulations is
sparse. There are a number of challenges in the 3D modeling, including the
computational complexity of the problem, a need for appropriate initial
conditions to ignite the inverse pinch, insufficient equation of state (EOS)
data for the working gas and the electrodes, and incorporating radiative
effects into the MHD model.

In a DPF, a working gas is charged, forming a plasma, that travels down an
anode surrounded by a cathode, and the speed that the plasma shock travels
down the anode is determined by the initial pressure and voltage in the
system. The Z-pinch occurs when the plasma shock gets to the end of the anode.
The energy available to the Z-pinch is maximized when the current through the
DPF is maximized, which means that the rundown time is one of the most
important quantities to properly simulate. If simulations correctly predict
rundown for a given initial pressure and voltage configuration, those settings
can be used in actual experimentation ensuring that the Z-pinch occurs with
maximum energy, which in turn results in maximum neutron yield. The neutron
yield can be determined experimentally using a Beryllium activation detector
[9] (for deuterium fusion), or a Praseodymium activation detector [10] (for
deuterium-tritium fusion). Further, when deuterium, or deuterium-tritium
mixtures are used as a fill gas, the neutron yield has a power-law
relationship [11] with the maximum current, which means that it is also
important to be able to faithfully simulate the maximum current. For these
reasons, most of our analysis centers on comparing experimentally-measured
current waveforms to the simulated current waveform.

In this work, we present simulation results of a fully 3D MHD model of a DPF
using the ALEGRA [12] multiphysics code developed at Sandia National
Laboratories. The simulations are run in both 2D and 3D, and the results are
compared to each other as well as to experiments run in the DPF lab at
National Security Technologies, LLC. The predictions of the 2D and 3D
computations are qualitatively and quantitatively different, as the 2D
simulations show systematically lower inductance, which results in
systematically higher currents but unrealistically fast rundown times. The 3D
simulations predict lower maximum current values but accurately represent the
true rundown time shown in experimental data. This demonstrates that, despite
the symmetric geometry of the machine, there are three-dimensional effects
present in the MHD physics that must be accounted for in order to faithfully
predict the outcome of DPF experiments.

##  II DPF Physics and the Experimental Setup

The DPF used in our experiments, and the geometry of which was modeled in the
simulations, is formed of coaxial electrodes in a rarified deuterium
atmosphere (about 7 Torr). A two-stage Marx capacitor bank is charged, and
when discharged, breaks down the gas, forming a shock and starting the
“inverse pinch” phase of the reaction, in which the gas expands outward from
the anode to the cathode bars (see Fig. 1). Once the gas touches the cathode,
the “run-down” phase begins, and the plasma moves up the anode until it
Z-pinches at the top of the anode. These two phases, the inverse pinch and
rundown, are approximately governed by magnetohydrodynamics and are the
components of the reaction that are studied and simulated.

###  II-A DPF Geometry and Setup

Fig. 1 shows the DPF setup used in our experiments. The outer electrode
(cathode) is formed of 24 copper bars, 0.375 inches thick and 30.75 inches
tall, in a ring with an inside diameter of 6 inches. The cathode is at ground
potential, and its bars are shorted at the top with a ring of copper. The
anode is a hollow copper tube with an outer diameter of 4 inches that stands
23.6 inches above the ground plane, capped with a hemisphere. The vacuum
chamber is 1 foot in diameter, and roughly 6 inches taller than the cathode
cage. A Pyrex insulator tube, which is about 0.5 inches thick and stands 8.63
inches above the cathode base, separates the anode and cathode.

Figure 1: A rendering of the DPF used for the models and experiments. The
anode is the dome-topped cylinder in the center; the cathode “cage” is the
collection of rectangular bars that surrounds the anode; and the insulator is
visible through the bars, at the bottom of the cathode cage. The vacuum
envelope is represented as the tall cylinder that surrounds the cathode and
anode, and the ground plate is the flat cylinder at the bottom. Some support
features, such as the cathode top support ring, have been left out of the
drawing.

The DPF is driven by a two-stage Marx capacitor bank, which is connected to
the plasma focus tube by 36 coaxial cables. The total capacitance of the bank
(when configured for discharge) is 432 $\mu$F, and the maximum total voltage
in discharge configuration is 70 kV, which makes the maximum stored energy of
the bank 1 MJ. The plasma shock is driven by an external circuit, and the
circuit model used in the MHD simulations is shown in Fig. 2. The discharge
switch in the circuit represents a collection of eight rail-gap switches that
are simultaneously triggered by a single spark gap. The series inductance
represents the transmission lines that feed power to the plasma focus tube and
was determined empirically by fitting exponentially-dampened sine waves to the
experimental data. The 10 nF capacitor in series with the small resistor
represents the imperfect capacitance of the terminal plates and the
transmission lines that supply the power to the plasma focus tube. The 120
$\Omega$ resistor in parallel with the plasma focus tube is the equivalent
parallel resistance of the safety resistors.

Figure 2: The equivalent DPF discharge circuit used in the MHD models. The top
connection of the DPF is to the anode, and the bottom connection is to the
cathode. The switch shown is a collection of eight triggered rail-gap
switches.

###  II-B Faraday Current Diagnostic

The discharge current diagnostic is important for comparing simulations to
experiment, because it allows for the measurement of both the maximum current
through and rundown time of the DPF. As was noted above, these are two of the
most important quantities for simulations in order to accurately predict
neutron yield and to ensure that the maximum current runs through the system
at the time of the Z-pinch. On the NSTec DPF, the discharge current is
measured with a Faraday rotator [13], which uses the magnetically-induced
linear polarization rotation in quartz fibers to measure the current in a
circuit. The fiber is wound in a circular fashion around the anode, an
orientation that causes the fiber to follow the direction of the magnetic
field, allowing it to accurately measure the current. In (1), the polarization
rotation angle, $\Theta$ (in radians), is related to the permeability of the
vacuum, $\mu_{0}$, the current, $I$, the Verdet constant of the fiber, $V$,
and the number of loops the fiber makes around the anode $n$. The interaction
of the magnetic field, $\vec{B}$, with an element of the fiber’s length,
$\vec{d\ell}$, is then integrated over the path of the fiber around the anode.
In our arrangement, the fiber wraps around the anode $n=5.25$ turns. This path
is denoted by $\xi$, and since the fiber is either parallel with the magnetic
field or perpendicular to it, the rotation angle is

$\Theta=V\int_{\xi}\vec{B}\cdot\vec{d\ell}=\mu_{0}nVI,$ (1)

in MKSA units. The Faraday rotator current diagnostic is perferred over other
discharge current measurements, such as the Rogowski coil, because the Faraday
rotator gives current measurements that are not dependent on calibration
factors, but rather on an easily measurable geometric quantitiy: the number of
turns around the current to be measured. The only other factor that must be
determined is the Verdet constant for the fiber, which can either be measured
or obtained from a datasheet on the fiber. When properly set up, the Faraday
rotator is a reliable diagnostic.

Figure 3: Shown is a comparison of the current profiles for thirty-seven DPF
shots, all at the same voltage and pressure (37.5 kV and 7.28 Torr,
respectively). This demonstrates the consistency of the current produced by
the machine, as well as representative Faraday rotator data.

Fig. 3 shows experimentally-measured results from the Faraday rotator detector
for 37 DPF shots initiated with the same voltage and pressure. As can be seen,
the profiles are all nearly identical, which demonstrates both the shot-to-
shot consistency of the DPF and the reliability of the Faraday diagnostic. The
placement of the Faraday loop is important for understanding the current that
it measures. This is easier to show than it is to describe, so this is
included as fig. 4.

Figure 4: Shown is a diagram of where the Faraday coil is placed on the Gemini
DPF. In the cutaway, the red (darker) area represents the parts electrically
connected to the anode, and the violet (lighter) areas are considered to be at
cathode potential. The Faraday loop is shown as a hoop that is under the
center center conductor wires, and above the ground plate.

##  III Modeling and Simulation of the DPF

The modeling for this project was performed with Sandia National Laboratories’
ALEGRA-MHD code. ALEGRA is a finite-element, multi-material, arbitrary
Lagrangian-Eulerian (ALE) shock hydrodynamic code designed for parallel
computing. It uses an operator-split edge-element formulation to simulate
resistive MHD in 2D and 3D high-deformation shocked media and pulsed power
systems.

ALEGRA provides fine control over how the simulation is performed through a
text file known as the “input deck.” The primary purpose of the input deck is
to define the simulation geometry, material composition, and physics to be
modeled, though it also allows a user to request output and provide runtime
controls. Part of defining the physics of a simulation is setting up the
equations of state for the gases, which was among the most challenging aspects
of the problem. Tabular equations of state provide the most effective means of
modeling the thermodynamic state of material in the simulations, which may be
in the solid, liquid, gaseous, or ionized state. Tabular EOS models were made
available through ALEGRA’s interface to Los Alamos National Laboratory’s
SESAME [14] data. Tabular representations of the Lee-More-Desjarlais (LMD)
electrical and thermal conductivity model [15] were also used here. The LMD
model combines empirical data with inferences from quantum molecular dynamics
modeling and density functional theory (QMD-DFT) to provide a conductivity
representation that spans the transition between conducting and insulating
conditions and has proven quite successful in this “warm dense matter” regime
[16, 17].

There are two approaches to defining the geometry on which to simulate the
reaction using ALGERA. The first is using ALEGRA’s built-in functionality, by
declaring pre-defined volumetric shapes (like spheres, prisms, pyramids, etc.)
and defining the material composition of each shape. For example, many DPF
machines have cylindrical copper cathode bars, and it is possible to define in
the input deck a cylinder made of copper. The code then generates a 2D or 3D
unstructured mesh that overlays the defined shapes, allowing for mesh elements
to intersect more than one user-defined shape. The second approach to defining
the geometry is by performing the computation on a body-fitted mesh generated
using an external meshing tool. The body-fitted mesh method is a more accurate
method of describing the material in the simulation volume, but is more time-
consuming to set up than the geometric method. ALEGRA’s built-in functionality
was used to define the geometry and material composition in all of the
simulations shown in what follows, and boundary conditions are specified on
subsets of nodes or faces within the simulation volume.

In any simulation, it is important that the geometry and the physics being
simulated reflect the important geometries and physics of the experiment, and
both raise several concerns in our simulations. For example, it is not
necessary to include a faithful model of the vacuum chamber in our
simulations, since the plasma does not interact with the top of the chamber
during the simulation. Evaluating the physics being modeled, it is important
to understand that the plasma in the DPF reaction is driven by an external
electrical network, and a great strength of the ALEGRA-MHD code is that it has
a sophisticated built-in circuit solver that can be used to couple electrical
energy from user-specified circuits into the simulation volume.

Near the end of the simulation, just prior to the Z-pinch, the MHD simulation
begins to become unphysical because of its inability to represent certain
phenomena, such as, the kinetic instabilities which raise the plasma
resistivity. The simulation may run past the point of Z-pinch without
crashing, but the time-evolution of the simulation would be unphysical. Once
the MHD simulation begins to approach the Z-pinch, it is possible to transfer
the model state information to a particle-in-cell model to accurately simulate
the Z-pinch, which has been demonstrated by researchers at LLNL[18] and
SNL[8].

Setting up the initial condition for the simulation can also be tricky, since
the DPF’s starting state happens when the gas in the chamber breaks down in
the vicinity of the insulating sheath and becomes a ring-shaped plasma shock.
The breakdown of the gas is not covered by MHD physics, so the gas near the
insulator has to be initialized in an artificial state that will quickly
transition to the plasma shock known experimentally to exist in the DPF. We
have found that setting up a thin layer of extremely hot ($\sim 10^{6}$ K, and
therefore conductive) gas on the surface of the insulator results in the
simulation initiating a plasma shock without causing observable artifacts in
the time evolution of the simulation. Our experience with this initial
condition is that the thin layer of hot gas should be about as thin as the
Pyrex insulator and should touch both the anode and the cathode. Using this
initiation of the plasma shock results in temperatures that match data from
particle-in-cell calculations [8]. The artificially hot gas layer should
stabilize its temperature near the shock temperature, $\sim 10^{4}K$, within a
few solver timesteps (typically, about 20 nanoseconds).

###  III-A Two-dimensional Modeling

One and two dimensional models of the DPF are the most common in the
literature, frequently coming in two generic types: empirical models and
finite-element MHD models. The primary strength of empirical models, such as
RADPFV5.5de [19], or Scat95 [20] is that they give results that are often very
close to experimental data. A significant drawback is that they require
existing data to fit the model to, and can therefore only be used to simulate
devices that are quite similar to the devices that one has experimental data
for already. The second class are finite-element codes such as MHRDR [21] and
Mach2 [22] that perform MHD modeling in one or two dimensions or in axis-
symmetric 3D. Fully 3D MHD codes are not new, but can be hard to obtain and
are usually more difficult to operate than 2D codes. Results of fully-3D DPF
modeling are not well represented in the literature.

In many experimental realizations of the DPF, the cathode is a cylinder
composed of metallic bars. While the plasma shock is traveling down the tube
during rundown, the plasma spills outside of the anode-cathode gap through the
spaces between the bars, a phenomenon that can be seen in the framing camera
picture in Fig. 5. This can be approximated in 2D simulations, but in practice
it is difficult to predict the time evolution of an actual DPF using these
approximations.

Figure 5: Shown is a framing camera picture of the plasma rising up the
electrodes in a DPF. The top of the anode can be seen as the dark disk in the
middle of the bars, which can be seen at the top. The lower region of the
chamber is bright due to the plasma, which has escaped the anode-cathode gap
and surrounded the bars.

Simulations in 2D also impose other geometric constraints: the current density
must either be completely in the radial-axial plane, or perpendicular to it,
but not both at the same time. This precludes modeling situations that may
have currents flowing helically, or situations in which the current may be
flowing asymmetrically or off-axis. These restrictions on 2D simulations are
often not of concern to investigators, who may want to ensure symmetry and
simplicity in their experiments in order to simplify the data they collect.
Nonetheless, the most general, physically realizable simulations are essential
for complete understanding.

The 2D simulation that was run in ALEGRA made the ad hoc assumption that there
was a lower density floor below which the material was assumed to have no
electrical or thermal conductivity. The floor was set at a density of
2.5$\times$10-4 kg/m3, which was was necessary in order to eliminate
unphysical behavior in the simulation. The LMD model for deuterium plasma,
shown in Figure 6, attributes moderate electrical conductivity to the plasma
at this density for temperatures higher than approximately 1 eV, and this
behavior is suppressed here, in order to cause the plasma shock to travel
properly down the anode of the DPF.

Figure 6: The Lee-More-Desjarlais (LMD) electrical conductivity model for
deuterium.

###  III-B Three-dimensional Modeling

Fully 3D magnetohydrodynamic codes are available, such as ALEGRA and NIMROD
[23], among others, and the benefit of using these codes is that the electrode
geometry in the DPF can be modeled and simulated. The ability to investigate
the effects of electrode structures without presupposing symmetries allows
investigators to gain deeper understanding into how the plasma shock evolves
over time. The primary drawback of 3D MHD modeling is the computational
complexity of solving the MHD equations on large meshes, requiring computer
clusters to perform simulations in a reasonable amount of time. Similar to the
2D modeling, the 3D modeling in ALEGRA required the imposition of a density
floor at 2.5$\times$10-4 kg/m3 in the electrical and thermal conductivity
models.

Fig. 7 shows a representation of the material density midway through the
rundown phase of the DPF in both the 2D and 3D simulations. The 2D modeling
assumes axial symmetry, and the symmetry axis in both simulation volumes is on
the far left-hand side. In the both simulations, the cathode bars are
represented as rectangles on the right of the simulation domain, and the
plasma shock travels from the bottom of the image to the top of the image,
where it Z-pinches slightly above the hemispherical anode top. The plasma
cannot flow around the cathode bars in the 2D simulation as it can in the 3D
simulation. The 3D simulation simulates the entire gas volume of the chamber,
and it can be seen in Fig. 7 that the plasma is slightly slower and less dense
in the 3D simulation than in the 2D simulation. The larger inductance results
in longer rundown times and lower maximum current as compared to the 2D
simulation, if all other system parameters are equal.

Figure 7: This is a comparison of the density on a slice through the
simulation volume at about 3 microseconds. The image on the left shows the 2D
simulation where the plasma cannot flow around the bars, and the image on the
right shows the 3D simulation, where the plasma can flow around the bars.
Since the plasma flows around the bars in the 3D simulation, it also affects
the external impedance as a function of time for the simulation volume.

##  IV Experimental and Predicted Results Comparison

Fig. 8 shows the current profiles from a 3D simulation, two 2D simulations,
and from the Faraday diagnostic of an actual experimental run. The “simulated
current” in 2D (red, dashed line) and 3D (yellow, dash-dotted line) were both
initiated using the same voltage and pressure as the experimental data. The 2D
simulation predicts a peak current of 2.08 MA, which differs from the peak
current measured by the Faraday diagnostic (2.17 MA) by only 4%. The 3D
simulation systematically predicts lower peak currents, due to the higher
inductance of the plasma escaping the cathode bars, and estimates peak current
at 1.82 MA, an error of approximately 16%. Thus, for estimating peak current,
the 2D simulation is more accurate than the corresponding 3D simulation.

The more important quantity of interest, however, is the duration time of the
rundown phase of the reaction. The end of the rundown phase is defined at the
point of maximum derivative of the current profile, which is 6.96 $\mu$sec for
the experimental data. The 3D simulation predicts 6.69 $\mu$sec, an error of
less than 4%, whereas the 2D simulation predicts a rundown time of 5.59
$\mu$sec, an error of almost 20%, showing that the 3D simulation vastly
outperforms the 2D simulation in predicting this quantity.

The two primary inputs into the ALEGRA simulation that we have discussed so
far are the initial voltage and pressure, since these are the adjustable
initial conditions of the DPF machine. The simulation allows for other
quantities to be tweaked as well, though, and it is natural to adjust the
model parameters to attempt to match the experimental data more closely. This
is difficult with the 3D simulations, since they are too computationally
intensive to “tweak” parameters one at a time and analyze the changes in
output. For the 2D simulations, however, this is possible, and Fig. 8 shows
the results of a 2D simulation (“tweaked current,” pink dotted line) that was
designed to match the rundown time of the experimental results. This required
that the series inductance be adjusted from 25 nH to 28.2 nH. Note that there
is no experimental justification for this change, it is done just to show that
the true rundown time can be achieved with a 2D simulation. This simulation
does not outperform the 3D, however, since the 2D simulation matching the
experimental rundown time results in a far inferior peak current measurement.
Thus a small improvement in rundown time over the 3D gives a large degradation
in peak current.

Figure 8: Shown is a comparison of the simulated 2D current (red, dashed line)
to the current measured by the Faraday probe (blue, solid line) and to the
simulated 3D current (yellow, dash-dotted line) and a Scat95 simulation
(green, dotted line). The 2D simulation underestimates peak current and
severely underestimates rundown time. The 3D simulation also underestimates
peak current but more faithfully predicts rundown time. Also shown is a 2D
simulation (pink, dotted line) whose input parameters are adjusted to give a
rundown time similar to the experimental data. In order that the 2D simulation
match the experimental rundown time, the peak current is severely
underestimated. The Scat95 simualtion shows better agreement, however,
requires iterative adjustment of parameters to already existing experimental
data.

Naturally one should use simulation input parameters that represent the
experiment being simulated as faithfully as possible, and the story of Fig. 8
is that 2D simulations of a DPF can give quite good results when the peak
current is the quantity of interest. When the rundown time is of interest,
which is more often the case, it is necessary to use the fully 3D simulation
to accurately predict the rundown time of an experiment. For comparison, a
Scat95 simulation was iteratively adjusted to match the experimental data.
While the agreement is good for this simulation, the parameters used in the
match are only good matches for geometries and setups that are close to this
particular case. Otherwise Scat95 achieves results that are similar to the 2D
MHD simulations.

##  V Conclusions

In this work we have presented results of fully 3D predictive simulations of a
Dense Plasma Focus, using the ALEGRA MHD code from Sandia National
Laboratories. As opposed to 2D or axis-symmetric 3D simulations, the fully 3D
models more faithfully predict the duration of the rundown phase of the DPF,
which is essential for ensuring that the maximum current runs through the
system at the time of Z-pinch, which is required to accurately predict neutron
yield. The 2D simulations are appropriate for predicting the peak current in
the DPF, but are not capable of matching both the peak current and rundown
time simultaneously.

## Acknowledgment

The authors would like to thank Aaron Luttman for helpful comments and
suggestions on the manuscript. We would also like to thank Chris Hagen for
providing support and encouragement for this project, as well as his insight
into the theoretical and experimental operation of the DPF.

Sandia National Laboratories is a multi-program laboratory managed and
operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin
Corporation, for the U.S. Department of Energy’s National Nuclear Security
Administration under contract DE-AC04-94AL85000.

This manuscript has been authored in part by National Security Technologies,
LLC, under Contract No. DE-AC52-06NA25946 with the U.S. Department of Energy
and supported by the Site-Directed Research and Development Program. The
United States Government retains and the publisher, by accepting the article
for publication, acknowledges that the United States Government retains a non-
exclusive, paid-up, irrevocable, world-wide license to publish or reproduce
the published form of this manuscript, or allow others to do so, for United
States Government purposes.

## References

  * [1] J. W. Mather, “Investigation of the High Energy Acceleration Mode in the Coaxial Gun,” Phys. Fluids. Suppl., vol. 7, S28, 1964.
  * [2] J. W. Mather, “Formation of a High-Density Deuterium Plasma Focus,” Phys. Fluids, vol. 8, no. 2, Feb., 1965, pp. 366-377.
  * [3] N. V. Filipov, T. I. Filipova, and V. P. Vinogradov, “Dense, High Temperature Plasma in a Noncylindrical Z-pinch Compression,” Nuc. Fusion, vol. 2, 577, 1962.
  * [4] V. A. Gribkov, L. Mahe, P. Lee, S. Lee, A. Srivastava, “Dense Plasma Focus radiation source for microlithography & micro-machining,” Proc. SPIE Microlithographic Techniques Integrated Circuit Fabrication II, vol. 4226, 2000, pp. 151-159.
  * [5] E. J. Lerner, S. K. Murali and A. Haboub, “Theory and Experimental Program for p-B11 Fusion with the Dense Plasma Focus,” J. Fusion Energy, vol. 30, no. 5, Jan., 2012, pp. 367-376.
  * [6] J. O. Pouzo, “Applications of the dense plasma focus to nuclear fusion and plasma astrophysics,” Plasma Science, IEEE Trans., vol. 31, no. 6, Dec., 2003, pp. 1237-1242.
  * [7] J. H. Gonzalez, F. R. Brollo, and A. Clausse, “Modeling of the Dynamic Plasma Pinch in Plasma Focus Discharges Based in Von Karman Approximations,” Plasma Science, IEEE Transactions on, vol. 37, no. 11, pp. 2178-2185, Nov., 2009.
  * [8] C. S. Kueny, D. G. Flicker, and D. V. Rose, “ALEGRA-HEDP Simulations of the Dense Plasma Focus,” Sandia National Laboratories, Albuquerque, NM, SAND2009-6373, 2009.
  * [9] T. J. Murphy, “A practical Beryllium activation detector for measuring DD neutron yield from ICF targets,” LA-UR-96-1649 Los Alamos Unclassified Report, 1996.
  * [10] B. T. Meehan, E. C. Hagen, C. L. Ruiz, G. W. Cooper, “Praseodymium activation detector for measuring bursts of 14 MeV neutrons,” Nucl. Instr. and Meth. A, vol. 620, 2010, pp. 397-400.
  * [11] O. Zucker, et. al., “The plasma focus as a large fluence neutron source,” Nucl. Instrm. and Meth., vol. 145, issue 1, Aug., 1977, pp. 185-190.
  * [12] W. J. Rider, A. C. Robinson, et al., “ALEGRA: An Arbitrary Lagrangian-Eulerian Multimaterial, Multiphysics Code,” Proc. 46th AIAA Aero. Sci. Meeting., Reno, NV, Jan., 2008.
  * [13] L.R. Veeser, G.W. Day, “Fiber-Optic, Faraday Rotation Current Sensor,” LA-UR-86-2084 Los Alamos Unclassified Report, 1984.
  * [14] S. P. Lyon, J. D. Johnson, “SESAME: The Los Alamos National Laboratory equation of state database.” LA-UR-92-3407 Los Alamos Unclassified Report, 1992.
  * [15] M. P. Desjarlais, “Practical Improvements to the Lee-More Conductivity Near the Metal-Insulator Transition,” Contrib. Plasm. Phys., vol. 41, 2001, pp. 267-270.
  * [16] M. K. Matzen, M. A. Sweeney, R. G. Adams, J. R. Asay, J. E. Bailey, et al., “Pulsed-power-driven high energy density physics and inertial confinement fusion research,” Phys. Plasmas, vol. 12, no. 055503, 2005.
  * [17] R. M. Lemke, M. D. Knudson, C. A. Hall, T. H. Haill, M. P. Desjarlais, J. R. Asay, “Characterization of magneticaly accelerated flyer plates,” Phys. Plasmas, vol. 10, no. 4, 2003, pp. 1092-1099.
  * [18] A. Schmidt, V. Tang, D. Welch, “Fully Kinetic Simulations of Dense Plasma Focus Z-Pinch Devices,” Phys. Rev. Letters, vol. 109, 205003 (2012).
  * [19] S. Lee, S. H. Saw, et al., “Characterizing Plasma Focus Devices - Role of the Static Inductance - Instability Phase Fitted by Anomalous Resistances,” J. Fusion Energy, vol. 30, no. 4, Aug. 2011pp. 277-282.
  * [20] R. Gribble, M. Yapuncich, W. Deninger, “SCAT95,” 2.0 Los Alamos National Laboratory, 1997.
  * [21] V. Makhin, B. Bauer, et al., “Numerical Modeling of a Magnetic Flux Compression Experiment,” Journal of Fusion Energy, vol. 26, 109-112, 2007
  * [22] M. H. Freese, “MACH2: A Two-Dimensional Magnetohydrodynamic Simulation Code for Complex Experimental Configurations,” AMRC-R-874, (1987).
  * [23] C. R. Sovinec, A. H. Glasser, et al., “Nonlinear Magnetohydrodynamics with High-order Finite Elements,” J. Comp. Phys., 195, 335 (2004).

B. T. Meehan received a B.S. in Physics from the United States Naval Academy
in 1995 and an M.S. in Applied Physics from Stanford University in 1997.
Currently he works with the Dense Plasma Focus Group at National Security
Technologies, LLC, modeling DPFs with magnetohydrodynamics codes.
---
J. H. J. Niederhaus holds a B.S. in Physics from the Virginia Military
Institute (2001), an M.S. in Nuclear Engineering from the Pennsylvania State
University (2003), and a Ph.D. in Engineering Physics from the University of
Wisconsin-Madison (2007). He is a computer scientist in the Computational
Shock and Multiphysics Department at Sandia National Laboratories.
---