# Statistical Learning for Best Practices in Tattoo Removal

Richard P. Yim1, Jamie Haddock1, Deanna Needell1 Correspondence:
<EMAIL_ADDRESS>
(Affiliations: 1Department of Mathematics, University of California, Los
Angeles)

###### Abstract

The causes behind complications in laser-assisted tattoo removal are currently
not well understood, and in the literature relating to tattoo removal the
emphasis on removal treatment is on removal technologies and tools, not best
parameters involved in the treatment process. Additionally, the very challenge
of determining best practices is difficult given the complexity of
interactions between factors that may correlate to these complications. In
this paper we apply a battery of classical statistical methods and techniques
to identify features that may be closely correlated to causes of complication
during the tattoo removal process, and report quantitative evidence for
potential best practices. We develop elementary statistical descriptions of
tattoo data collected by the largest gang rehabilitation and reentry
organization in the world, Homeboy Industries; perform parametric and
nonparametric tests of significance; and finally, produce a statistical model
explaining treatment parameter interactions, as well as develop a ranking
system for treatment parameters utilizing bootstrapping and gradient boosting.

Keywords:  Laser Assisted Tattoo Removal, Parametric Tests, Non-Parametric
Tests, Binary Classification, Logistic Regression, Gradient Boosting,
Bootstrapping

##  1 Introduction

Current best treatment practices for laser-assisted tattoo removal are ad hoc
at best, and the interactions between parameters involved in the tattoo
removal treatment process are not well understood. For instance, it is
believed that patient demographic factors may potentially affect treatment
outcomes (e.g., factors such as patient age and gender), and that these
factors are a source of variation in physiological processes between
individuals (i.e., skin reactivity and healing time vary demographically).
Additionally, much of the research surrounding tattoo removal is focused on
the clinical trial setting with an emphasis on discovering new innovations for
laser-assisted tattoo removal [1]. While there is considerable debate on laser
pulse duration standards, it still remains that there is no precise and
standard system for best practices that has been shown to be both optimally
effective and safe [13]. A reason for the lack of understanding of safe
treatment practices is a general lack of records of treament parameters
causing complications [9].

For our study we explore data recorded by Homeboy Industries, the largest gang
rehabilitation program in the world, which offers free tattoo removal
treatment services for former gang members seeking to remove gang-related
tattoos. We utilize various statistical methods and demonstrate some results
of inference relating to laser-assisted tattoo removal treatment procedures
performed in a practical clinical setting at Homeboy Industries in Los
Angeles.

In our data there are four types of complications that are associated with
tattoo removals: hyperpigmentation and hypopigmentation, or increased and
decreased skin pigmentation, respectively; scarring, visible tissue regrowth;
and keloids, scarring with excessive skin overgrowth in the treated area. We
emphasize that while the complication rate among the studied sample is low,
the main concern for all tattoo-removal practitioners is a commitment to doing
as little harm as possible. For this reason, our primary concern with the
available data is to simply understand complication occurrence and
nonoccurence in general, grouping all instances of complications together
without distinction. The main significance of this design is to understand the
tattoo removal process as one that may potentially discourage patients seeking
removal from continuing the treatment process, and cause unnecessary grief in
addition to any preexisting tattoo regret, as is commonly surveyed amongst
individuals with tattoos [16].

Generally, in this paper we apply a full range of inference techniques and
machine learning methods to develop an understanding of best tattoo removal
treatment parameters. The data measures treatment settings used on a tattoo
over multiple treatment appointments. Data is provided at both the patient
level and the treatment level for a given tattoo. With regards to methods
used, we begin with analyzing the significance of patient demographic factors
that are understood to be correlated with complication occurrence by utilizing
simple statistical tests. We then perform parametric and nonparametric tests
of significance for analysis of variance (e.g., Kruskal-Wallis, Wilcoxon rank
sum and randomization) and use Tukey’s honestly significant difference (HSD)
test to identify important sources of variation. Finally, we fit a logistic
regression model to identify statistically significant treatment parameters,
and perform multiple gradient boosting using decision trees to form a ranking
of treatment parameters that a health care practitioner should be most aware
of during laser-assisted tattoo removal procedures.

This paper is organized as follows: Section 2 details the dataset and provides
descriptions of the methods utilized; Section 3 presents results of the
methods employed; we conclude the paper in Section 4 with last remarks.

##  2 Data and Methods

The individuals that come through Homeboy industries are often former gang
members, and many of them have tattoos that may be offensive or attract
negative attention. In this section, we discuss the various datasets recorded
by Homeboy Industries. The final cleaned datasets that were used in our
analysis detailed tattoo level data and patient demographic data.

We also utilized many different statistical tests and methods, both parametric
and nonparametric. We detail the properties of these tests and their
requirements for legitimate application to the data. All computations were
done with R-4.0.3 [18]. We used the caret library as the primary driver for
applying the machine learning algorithms in our study. All computing was done
on a Linux machine running Ubuntu 20.04 with an Intel 10700K processor.

###  2.1 Patient Demographic Data

Including observations with missing values, there are a total of 2,118 tattoo
observations among 502 patients recorded in the data. We briefly detail the
various factors that were recorded at the patient demographic level in Table
1. We note the importance of the Fitzpatrick score—a measure developed in 1975
to classify skin tones [6]—as it is currently believed in the literature that
the efficacy of selected treatment parameters for laser-assisted tattoo
removal is dependent on the Fitzpatrick score.

Table 1: Table detailing patient demographic and patient-level factors of interest. Variable |  Description
---|---
Patient Age |  Age of the patient as of 10 June 2020 (integer-valued)
Sex |  Male or female (one-hot encoded as male/not-male, binary)
Ethnicity |  Hispanic/Latino, or not (one-hot encoded, binary)
Race |  Pacific Islander, American/Alaskan Indian, Black, Asian, Latino/Hispanic, white, multiracial, other (nominal)
Treatment Total |  Total number of treatment visits at Homeboy for tattoo removal (integer-valued)
Total Tattoos |  Total number of tattoos listed by patient (integer-valued)
Fitzpatrick Score |  Fairest tone (I) to deeply pigmented (VI) (ordinal)
Complications |  Indication of whether a patient ever experienced any complication (one-hot encoded, binary)

###  2.2 Tattoo Level Data

At the tattoo level, the data provided by Homeboy Industries includes
characteristics of tattoos that went through laser assisted tattoo removal. At
the tattoo level we perform statistical analysis on various factors grouped by
whether the observed tattoo ever experienced a complication in the treatment
sequence. We briefly detail the tattoo level features in Table 2. The data at
the tattoo level are important because it provides insight into the
significance of tattoo composition in whether a complication is likely to
develop as a result of undergoing laser treatment as recorded, as well as
challenge some assumptions of existing best practices.

Table 2: Table detailing tattoo level factors of interest. Variable |  Description
---|---
Category |  Tattoo location on body (nominal)
Age |  Age of tattoo (integer-valued)
Colors |  Whether the tattoo was black/blue or had other colors (one-hot encoded as black-blue/not-black-blue, binary)
Professional |  Whether the tattoo was done professionally, or by an amateur (one-hot encoded as professional/not-professional, binary)
Treatment Total |  Total number of treatments recorded for a tattoo (integer-valued)
Fitzpatrick Score |  Fitzpatrick skin tone rate: fairest tone at I to deeply pigmented at VI (ordinal)
Complications |  Indication of whether the tattoo ever experienced any complication at all (one-hot encoded, binary)

###  2.3 Treatment Level Data

For each tattoo at the treatment level we detail four main parameters involved
in laser-assisted tattoo removal, parameters that are actual tattoo-removal
laser settings: fluence, spot size, wavelength and frequency. We briefly
detail the treatment-level features in Table 3. Apart from observing
complication occurrence/non-occurrence as a response variable in our analysis
with respect to tattoo-level factors and patient demographic features, we also
study the variation of treatment parameters (e.g., fluence, spot size)
selected by practitioners given a particular tattoo characteristic such as
tattoo color and age. The intention for viewing treatment parameters as a
response characteristic is to gauge whether current clinician practices are at
all particular to tattoo-level characteristics at a statistically significant
level.

In our study, some additional design choices were made with regards to
studying the treatment parameters over a given treatment sequence. Since the
composition of tattoos have great variation between and within different
tattoo characteristics (e.g., color, ink composition and size) the number of
treatments recorded varied greatly between tattoos. For example, a treatment
sequence may have lasted over 10 weeks for one recorded tattoo, another may
have had only a single appointment. Thus, two variations of variable
transformation on the laser treatment parameters were made.

Table 3: Table detailing treatment level variables of interest. Variable |  Description
---|---
Fluence |  Laser heat and energy intensity measured in joules/cm2 (continuous)
Spot Size |  Laser spot radius measured in millimeters (continuous)
Wavelength |  Laser wavelength at two levels (532nm and 1064nm, binary)
Frequency |  Laser frequency measured in Hertz at two levels (5Hz and 10Hz, binary)
Treatment Day |  Days since first treatment (integer-valued)

The first variation is to simply record the mean and standard deviation of the
above treatment settings for the entire period in which a tattoo had undergone
treatment; this variation was used to study current clinician practices at
Homeboy Industries to see how responses of laser settings changed (i.e., how
clinicians selected different treatment parameters) on average as a result of
varying tattoo characteristics. The second variation also records the means
and standard deviations of the distribution of treatment parameters in
sequence, but only up to the first instance of a complication occurring;
naturally, if a complication never occurs in a tattoo, the mean and standard
deviation of applied laser settings are computed over the entire treatment
period. Figure 1 shows time series of fluence (in red), one of the four
possible laser parameters, over different tattoos with data from the full time
series up until the first-arrival of a complication—as well as change the of
laser fluence between treatments (in blue), to be detailed below. Note the
irregular duration of time intervals between treatments, as well as clear
nonlinear fluctuations of fluence chosen by clinicians.

Figure 1: We show fluence levels fluctuating by clinician application of
laser-assisted tattoo removal. The plots in the first row are of fluence
settings for a tattoo that does not develop any complications, while the plots
in the second row are of one that has experienced complications. The red time
series graphs represents the recorded fluence levels, while the blue time
series graphs represent the corresponding forward differenced time series of
fluence levels.

Furthermore, as an extension of the first complication arrival data, we create
additional variables of means and standard deviations with forward finite
differences. For a given variable, representing a time-series treatment
parameter, we compute mean and standard deviations for

$\Delta x_{k}=x_{k}-x_{k-1},$

where the terminating point in the time series, $x_{t}$, is dependent on
whether we are observing first-arrivals of complications or the entire
treatment sequence for a treated tattoo. For example, if we are observing
laser fluence for first-arrivals data we produce a forward finite difference
variable $\Delta x$, where the time series terminates at the index in the
sequence corresponding to the treatment where a complication was first
reported. Consequently, when computing the mean and standard deviations with
respect to this forward finite difference our expressions simplify as

$\overline{\Delta x}=\frac{1}{N_{t}-1}\left(\sum_{k=2}^{N_{t}}\Delta
x_{k}\right)=\frac{1}{N_{t}-1}\sum_{k=2}^{N_{t}}\left(x_{k}-x_{k-1}\right)=\frac{1}{N_{t}-1}\left(x_{N_{t}}-x_{1}\right),$

where we are essentially measuring the average likelihood of a given treatment
parameter setting from the first appointment to the terminal appointment
(i.e., the full treatment sequence for original treatment-level dataset in its
entirety, and the truncated treatment sequence for the first-arrivals
dataset). Similarly, for standard deviations of the differenced time series:

$\hat{\sigma}_{\Delta x}=\sqrt{\text{Var}{\left(\Delta
x\right)}}=\frac{1}{N_{t}-1}\sum_{k=2}^{N_{t}}\left(\Delta
x_{k}-\overline{\Delta
x}\right)^{2}=\frac{1}{N_{t}-1}\left(\sum_{k=2}^{N_{t}}(\Delta
x_{k})^{2}\right)-\left(\overline{\Delta x}\right)^{2},$

where $N_{t}$ is the number of treatments given to a particular tattoo. The
purpose of this variable differencing is to eliminate any apparent trend; a
first order forward differencing effectively applies a linear filter to our
time series [17]. The plots of the time series in Figure 1 demonstrate the
variation of $x_{k}$ and $\Delta x_{k}$ with respect to laser fluence. Tattoo
removal clinicians generally aim to gradually increase the fluence over the
treatment sequence, the effect of this difference is to essentially flatten
this linearly trending process, and characterize the distribution of applied
laser treatment settings into statistical parameters (e.g., mean and standard
deviation).

###  2.4 Statistical Inference

We applied both parametric and nonparametric statistical tests of significance
to identify statistically significant factors correlated to complication
occurrence in our data (in all our tests we use $\alpha=0.1$). We were
interested in observing differences in proportions across different patient
demographics and tattoo level factors. We applied the $z$-test for two sample
proportion testing as well as the chi-square test of independence for more
than two sample proportions. The hypothesis statement generally for $k$-sample
proportions is

$\displaystyle H_{0}:p_{a}-p_{b}=0$ $\displaystyle H_{a}:p_{a}-p_{b}\neq 0$

where $p_{a}$ and $p_{b}$ represent distinct sample proportions, pairwise,
across $k$-samples. The necessary assumptions for these tests were found to be
met before their application: independent samples and observations, sufficient
sample sizes, and mutually exclusive categories in variables [15].

Along with these statistical descriptions of proportions, we additionally were
interested in observing current distribution differences in removal practices
and laser treatment settings. We use Welch’s two-sample t-test to observe
parametric differences of averages between different treatment settings by
whether a tattoo experienced a complication or not [20, 21]. When the
necessary assumptions for a t-test failed—observation independence, sufficient
sample size and normality—we resorted to using the Wilcoxon rank sum test and
the randomization test. We strictly applied the Wilcoxon rank sum test to
responses that were categorical-ordinal on independent observations [14]. The
null hypothesis for the Wilcoxon rank sum test is that the difference in
central tendency of the distributions in our data is randomly occurring, with
the alternative being that difference in central tendency is not random, at a
statistically significant level.

In our implementation of the randomization test we were interested to see if
the measured difference of sample averages for settings applied to tattoos
that developed complications versus those that did not were of random chance.
The algorithm for the randomization test is as follows:

  1. 1.

Compute and store the true sample average difference.

  2. 2.

Randomization process:

    1. (a)

Remove labels of sample designation corresponding to each value; randomly
permute data. This step corresponds to null hypothesis of observation classes
as random.

    2. (b)

Relabel observations, partitioning data proportionally into new classes, and
measuring the differences of the average of the permuted data.

    3. (c)

Store value as part of an empirical distribution of randomized sample
statistics.

    4. (d)

Iterate for a sufficiently large $N$.

  3. 3.

Compute the empirical proportion from the randomization procedure of
differences of averages that are greater than or less than the true difference
(two-tailed variation) for an empirical $p$-value.

The null hypothesis of the randomization test is that the statistic of
interest is produced by random chance, with the alternative being that the
observed sample statistic is not randomly produced at a statistically
significant level. Both the Wilcoxon rank sum and randomization tests allow
for an intepretation of differences between factors that are correlated to
complication rates, albeit with less statistical power being nonparametric
with relative sample sizes between complication and noncomplication tattoos
being heavily unbalanced, 2,000 to 118, respectively [7].

Finally, in observing any true variation across factors with multiple levels
we applied one-way analysis of variance (ANOVA) where parametric assumptions
for normal probablity distributions were satisfied, and a Kruskal-Wallis test
for one-way ANOVA where these assumptions were violated. We verified the
assumptions of normality, linearity and homoscedasticity [2, 4]. The
hypotheses for one-way ANOVA are

$\displaystyle H_{0}$ $\displaystyle:\mu_{1}=\mu_{2}=\cdots=\mu_{k}$
$\displaystyle H_{a}$ $\displaystyle:\mu_{i}\neq\mu_{j},\;\text{for
some}\;i\neq j$

where the null is that the difference of averages between different categories
are not statistically significant, with the alternative being that there
exists a pair samples ($i$ and $j$) in the data where a statistically
significant source of variation in the data is present. In order to locate the
source of variation where shown to be significant, we applied Tukey’s HSD test
as a multiple comparison among the multi-level factors for which one-way ANOVA
was applied [19]. In particular, we used the Tukey-Kramer HSD test to account
for unequal samples when computing the studentized range distribution
statistic for the corresponding confidence interval

$\bar{y}_{i\cdot}-\bar{y}_{j\cdot}\pm\frac{q_{\alpha;k;N-k}}{\sqrt{2}}\hat{\sigma}_{\epsilon}\sqrt{\frac{1}{n_{i}}+\frac{1}{n_{j}}}\frac{}{}$

where $i$ and $j$ represent distinct sample distributions from different
populations; $q_{\alpha;k;N-k}$ follows a studentized range distribution with
parameters $\alpha$, representing the type II error; $k$ is the number of
populations or levels in a factor; and $N-k$ is degrees of freedom with $N$
being total number of observations; and without loss of generality,
$\bar{y}_{j\cdot}$ and $n_{j}$ represents the sample average for the feature
of interest and sample size, respectively [11]. Additionally,
$\hat{\sigma}_{\epsilon}$ denotes the square root of the overall variation of
the feature of interest in our total population.

When our parametric conditions failed (e.g., sample distribution normality and
homoscedasticity), we resorted to applying the Kruskal-Wallis test for one-way
nonparametric ANOVA [12]. The null hypothesis of this test is that the medians
of our populations are equal, and the alternative is that the there is at
least one source of random variation among our populations with a median that
is statistically significant among some pair of populations (e.g., median
fluence of settings applied to tattoos on the face versus tattoos on the neck
or those on the upper extremeties).

###  2.5 Statistical Learning

In addition to performing inference, we were interested in developing
statistical models of our data as related to complication occurrence (and
nonoccurrence) across the first-arrivals tattoo-level dataset. In our data for
complication occurrence, the actual response corresponds to a binary
classification problem where complication occurrence in a given tattoo
treatment sequence corresponds to a “true” outcome and complication
nonoccurrence corresponds to a “false” outcome. Generally, logistic regression
models a probabilistic relationship of the occurrence of some response along
with its predictors [10]. More precisely, logistic regression models the
following:

$p(X_{1},X_{2},\ldots,X_{n})=P(Y=1|X_{1},X_{2},\ldots,X_{n}),$

where in our case we have a binary label response of complication one-hot
encoded as “1” and complication nonoccurrence encoded as “0”; and there are
$n=10$ predictors, after heuristically removing highly correlated variables
using an absolute correlation threshold of 0.6. The $p$ represents the
probability that a given observation with its corresponding inputs are of an
outcome relating to tattoo complication. The precise formulation of the
logistic regression model is

$p(X_{1},X_{2},\ldots,X_{n})=\frac{e^{\beta_{0}+\beta_{1}X_{1}+\cdots+\beta_{n}X_{n}}}{1+e^{\beta_{0}+\beta_{1}X_{1}+\cdots+\beta_{n}X_{n}}}$

The coefficients of the logistic regression model are estimated using
maximimum likelihood estimation (MLE) of the log-likelihood function [3].
These coefficients asymptotically follow a normal distribution, from which we
can infer statistical significance and evident dependence on the corresponding
variables to the response in the data. We use this process to identify
statistically significant predictors of treatment parameters corresponding to
complications (again, with an $\alpha=0.1$).

In addition to modelling our data to find statistically significant factors
corresponding to complication occurrence, we wanted to produce a ranking of
statistically significant variables in our data. We do this using decision
trees, gradient boosting and bootstrapping [8, 5]. Decision trees are machine
learning models that recursively stratify (split) a feature space—tattoo
treatment parameters in our case—into regions that are estimated to be
numerically associated to a given response variable. Boosting is a general
statistical/machine learning paradigm that involves bootstrapping a datatset
into many copies, and sequentially fitting multiple machine learning models to
the dataset. In this boosting procedure on decision trees there are three
tuning parameters: $B$, the number of trees; $\lambda$, a shrinkage parameter
controlling the learning rate of the boosting method; and $d$, the number of
splits in a decision tree, or variables to consider in the feature space
stratification process.

Boosting with decision trees on a given dataset can produce ranks of variable
importance through the mean decrease in Gini index, a measure of total
variance across $K$ classes, formulated as

$G=\sum_{k=1}^{K}\hat{p}_{mk}(1-\hat{p}_{mk}),$

where $\hat{p}_{mk}$ represents the proportion of training observations in the
$m$th stratified region corresponding to the labels of class $k$. The rank of
variable importance is measured by the mean decrease in the Gini index over
all the decision trees fitted in the boosting procedure, where the higher the
mean decrease in the Gini index for that variable, the more important it is.
From this procedure we additionally bootstrap our dataset by fitting these
gradient boosted trees in our dataset 300 times, with each fit utilizing tuned
hyperparameters from repeated $k$-fold cross validation (8 folds for 2
repeats). We then produce an empirical distribution of ranks of variable
importance for features consdiered in our model; the most frequent rank of a
particular feature is then assigned as the feature’s empirical importance
rank. The final rankings of these features then decide which features tattoo
removal practitioners should be most aware of in the tattoo removal treatment
process.

##  3 Results

We tabulate our results of statistical inference as well as present the
results of our machine learning models. Throughout we discuss results and the
particular application of our methods to data.

###  3.1 Statistical Descriptions

Below are tables corresponding to two-sample (Table 4) and $k$-sample
proportion tests (subsection 3.1). We measure complication rates across
different factors at the patient demographic and tattoo levels. In our tables,
“Tattoo Complication” refers to whether a single tattoo ever experiences a
complication (i.e., complication occurrence/non-occurrence) within a treatment
sequence spanning multiple days.

Table 4: Details of two-sample proportions and p-values for $z$-test. Statistically significant factors determining significant proportion differences in complication rates differences are marked by the asterisk. The “True” column represents proportions that are characterized by the factor of interest; “False” column represents proportions that do not have the given factor characteristic. Response Proportion by Factor | True | False |  $p$-value
---|---|---|---
Treatment Completion by Complication | 0.14407 | 0.07300 | 0.0083*
Tattoo Complication by Colored Tattoo | 0.10857 | 0.04996 | 0.0035*
Tattoo Complication by Professional | 0.09351 | 0.04291 | 0.0010*
Complication by Sex (Male/Female) | 0.09964 | 0.11764 | 0.6163
Complication by Patient Median Age | 0.13061 | 0.08560 | 0.1381
Complication by Patient Ethnicity | 0.12010 | 0.06723 | 0.1452
Complication by Tattoo Median Age | 0.06379 | 0.08560 | 0.5562
Complication by Patient Fitzpatrick Score | 0.12000 | 0.10287 | 0.5562
Table 5: Details of $k$-sample proportions and p-values for chi-squared tests
of independence. Statistically significant factors determining significant
proportion differences in complication rates are marked by the asterisk.

Reponse by Factor | Parameter | 1st Quartile | 2nd Quartile | 3rd Quartile | 4th Quartile
---|---|---|---|---|---
| Quartile | 4 | 7 | 12 | 65
| $n_{i}$ | 122 | 116 | 117 | 122
Complication Rate by Total # of Treatments Tattoo ($p$=6.51e-10*) | $\hat{p}_{i}$ | 0.02417 | 0.03349 | 0.05181 | 0.12277
| Quartile | 4 | 7 | 12 | 100
Complication Rate by Total # of Tattoos per Patient ($p$=0.08395*) | $n_{i}$ | 92 | 70 | 87 | 100
| $\hat{p}_{i}$ | 0.11957 | 0.18571 | 0.08046 | 0.07000
| Quartile | 11 | 24 | 52 | 72
| $n_{i}$ | 122 | 116 | 117 | 122
Complication Rate by Total # of Treatments Patient ($p$=1.35e-06*) | $\hat{p}_{i}$ | 0.04918 | 0.06034 | 0.09402 | 0.24590
| Quartile | 1.4375 | 1.9134 | 2.5445 | 4.5773
Complication Rate by Mean Fluence (J/cm2) ($p$=3.532e-05*) | $n_{i}$ | 508 | 505 | 509 | 509
| $\hat{p}_{i}$ | 0.01771 | 0.04554 | 0.07466 | 0.07662
| Quartile | 3 yrs | 7 yrs | 12 yrs | 40 yrs
| $n_{i}$ | 292 | 306 | 233 | 253
Tattoo Complication Rate by Tattoo Age ($p$=0.1316) | $\hat{p}_{i}$ | 0.03425 | 0.071895 | 0.05150 | 0.075099
| Quartile | 29 yrs | 34 yrs | 41 yrs | 80 yrs
| $n_{i}$ | 148 | 109 | 128 | 117
Patient Complication Rate by Patient Age ($p$=0.2807) | $\hat{p}_{i}$ | 0.09459 | 0.07339 | 0.14844 | 0.11111

In our application of Welch’s $t$-test, Wilcoxon’s rank sum and randomization
tests, we observe the distribution differences of settings applied to the
tattoos over treatments up until the first arrival of a complication in the
treatment sequence; if a tattoo never experiences a complication the
statistics of treatment parameters is computed across all treatments for that
tattoo. Table 6 shows the results from our t-test as applied to distributions
of time-series statistics (mean and standard deviation) computed over our
settings that satisfies some heuristic observations of normality; Table 7
shows the results of Wilcoxon rank sum and randomization tests over
distributions of sample statistics that did not satisfy normality and follower
highly irregular distribution shapes (i.e., multimodal and skewed).

Table 6: We present the sample averages of these settings along with the corresponding $p$-values from Welch’s $t$-test. We find no statistically significant treatment settings among these parameters. Treatment Parameter | Complication | No Complication |  $p$-value
---|---|---|---
Mean Fluence | 1.958739 | 2.013944 | 0.5713
Mean Differenced Fluence | 0.1848074 | 0.1729800 | 0.8181
Mean Spot Size | 4.875772 | 4.942978 | 0.4402
Mean Differenced Spot Size | -0.06056931 | -0.08294852 | 0.6542

Table 7: We observe nonparameteric differences of statistics of treatment settings that are irregularly distributed with respect to the “first-arrivals” of complications data. For the Wilcoxon rank sum test we compute the sample median, and for the randomization test we compute the sample average. For each treatment parameter we present the results in two rows for sample median ($\hat{\eta}$, Wilcoxon rank sum) and sample mean ($\hat{\mu}$, randomization). Treatment Parameter | Statistic | Complication | No Complication |  $p$-value
---|---|---|---|---
| $\hat{\eta}$ | 1064 | 1064 | 9.197e-06*
Mean Wavelength | $\hat{\mu}$ | 1010.425 | 1050.169 | 7e-04*
| $\hat{\eta}$ | 0 | 0 | 0.05684*
Mean Differenced Frequency | $\hat{\mu}$ | 0.23336582 | 0.03502956 | 1
| $\hat{\eta}$ | 0.80475 | 0.72648 | 0.05307*
Spot Size Standard Deviation | $\hat{\mu}$ | 0.8165863 | 0.6899222 | 1
| $\hat{\eta}$ | 0 | 0 | 1.784e-06*
Mean Differenced Wavelength | $\hat{\mu}$ | 4.663210 | 1.571152 | 1
| $\hat{\eta}$ | 9.375 | 9.375 | 0.5299
Mean Frequency | $\hat{\mu}$ | 9.056429 | 9.145180 | 0.5633
| $\hat{\eta}$ | 61 | 71.5 | 0.5928
Average Days Between Treatments | $\hat{\mu}$ | 95.78247 | 102.58342 | 0.6128

Additionally, we have our results for current clinician practices and
approaches. Table 8 presents our results from the Kruskal-Wallis one-way ANOVA
for treatment parameters with distributions that violated the assumptions
required for parametric one-way ANOVA. We applied non-parametric one-way ANOVA
to study current variations of treatment parameters that were are found to be
of practical interest in the literature. In particular we were interested in
how clinicians applied different settings based on the tattoo age quartile.

Table 8: Sample averages by tattoo age quartile and Kruskal-Wallis $p$-values. All parameters listed were found to be statistically signficant (i.e., clinicians exercised distinct application of laser based on tattoo age). Quartile |  1st(3yrs) |  2nd (7yrs) |  3rd (12yrs) |  4th (40yrs) |
---|---|---|---|---|---
$n_{i}$ | 292 | 306 | 233 | 253 |  $p$-value
Mean Fluence ($\hat{\mu}_{i}$) | 1.96327 | 2.031156 | 2.019911 | 2.194109 | 0.001946*
Mean Spot Size ($\hat{\mu}_{i}$) | 4.888663 | 4.849143 | 4.857013 | 4.884559 | 7.314e-05*
SD Differ. Frequency ($\hat{\mu}_{i}$) | 2.374471 | 1.948435 | 0.08045977 | 1.948582 | 2.294e-05*

In our application of parametric one-way ANOVA for studying means, we found
particular interest in the application of average fluence with regards to the
location of the tattoo on the body. There indeed exists some statistically
significant variation when with how a clinician chooses a particular fluence
level with respect to tattoo location. Across tattoo locations on the body,
the distributions of fluence were found to be relatively constant in variation
and normally distributed. We obtained a $p$-value of 3.61-e05 for our one-way
ANOVA. Then, applying Tukey’s HSD test among tattoo location, we found that
tattoos applied to the upper extremities and face had the most significant
variation in average fluence applied—with sample averages of 2.09J/cm2 and
1.996J/cm2, respectively—among other tattoo locations on the body such as
lower extremities, back, chest, neck, head and abdomen.

###  3.2 Significant Factors

We report our machine learning results in this section. Table 9 shows the
numeric results of rank from our bootstrapped gradient boosting decision trees
as well as the coefficient estimates and corresponding $p$-values from our
logistic regression model. It should be noted that for both the logistic
regression and gradient-boosted decision tree models, we removed a few
variables from consideration. The removed variables were found to have high
absolute correlation ($>0.6$) with other variables in the first-arrival sample
statistics dataset of treatment parameters. The intention of reducing the
number of predictors was to simplify our model and reduce any potential
inflation of explained variation in our results.

With regards to the variable importance ranking system, with over 300
simulations of fitting gradient boosted trees, we produced a discrete
distribution of ranks for our variables, as detailed in Section 2.5. In cases
where the mode of ranks were tied between two factors, these ties were broken
by observing which rank mode was greater. For example, in Table 9, we found
that mean wavelength and mean differenced wavelength had rank modes of 1;
however, since the frequency of the rank mode for mean wavelength (100) was
greater than the rank mode of the mean differenced wavelength (90), we
reassign “total ranks” for mean wavelength as greater than the mean
differenced wavelength. Figure 2 displays a few rank distributions of variable
importance from our simulations.

Table 9: Results from boosted decision tree simulations, and coefficient estimates and $p$-values for logistic regression coefficients. Statistically significant logistic regression coefficients are asterisked. For our boosting procedure we show the ranking of important features corresponding to the mode of the variables’ frequency distribution. Ties in mode frequency rank were broken by comparing relative frequencies. Statistically significant logistic regression coefficients are asterisked. Treatment Parameter | Total Rank | Rank Mode | Estimate |  $p$-value
---|---|---|---|---
Mean Wavelength | 1 | 108 (1) | -0.004299 | 0.00118*
Mean Differenced Fluence | 2 | 84 (1 $\to$ 2) | 0.550567 | 0.33951
Mean Differenced Wavelength | 3 | 80 (1 $\to$ 2 $\to$ 3) | 0.004979 | 0.35263
SD Spot Size | 4 | 117 (4) | 1.294140 | 0.00609*
Mean Days Between Appts. | 5 | 93 (5) | -0.001966 | 0.35287
Mean Differenced Spot Size | 6 | 76 (5 $\to$ 6) | 0.449331 | 0.38701
Mean Differenced Frequency | 7 | 64 (6 $\to$ 7) | 0.466965 | 0.04217*
Mean Spot Size | 8 | 75 (7 $\to$ 8) | -0.154691 | 0.71596
Mean Fluence | 9 | 106 (9) | -0.137015 | 0.59433
Mean Frequency | 10 | 197 (10) | 0.014290 | 0.93119
Figure 2: Rank distributions of three variables with bars corresponding to
“total rank” hatched. We assign the rankings for each treatment parameter
based on the mode of their rank. We see that Mean Wavelength (Rank 1), Mean
Days Between Appts. (Rank 5), and Mean Fluence (Rank 9) are of descending
order of importance based on the location of the modes of ranks for each
variable from 300 simulations.

##  4 Conclusion

For our tattoo data we applied a whole range of methods from classical
statistical parametric tests to modern machine learning algorithms. We
characterized the time series treatment data across all parameters using
sample means and standard deviations, as well as having performed forward
finite differencing to characterize the variation in laser-assisted tattoo
removal treatment parameters between appointments. From the transformations
made on our variables, we applied both parametric and nonparametric tests of
significance to identify tattoo-level and patient demographic characteristics
strongly correlated to complication occcurence.

Below we list response-factor relationships that were found to be
statistically significant along with a brief details with regards to practical
implications:

  * •

Complication occurrence is related to decreased likelihood of tattoo removal
treatment completion. Patients may become ambivalent about continuing tattoo
treatments if a complication is experienced.

  * •

Tattoo removal complications are statistically more likely for tattoos done by
a professional tattoo artist as opposed to an amateur artist. This may be due
to the depth of the tattoo placement and/or ink utilized in the professional
tattoo process.

  * •

Tattoo removal complications are statistically more likely for tattoos with
color (e.g., red, green, yellow) as opposed to black/blue tattoos. This may be
due to the skin reaction sensitivity in tattoos with certain pigments.

  * •

Greater total number of treatments on a tattoo, number of tattoos on a
patient, and number of treatments a patient has undergone are all related to
increased complication rates. With regards to tattoo age, clinicians apply
statistically greater average fluences and have greater variability in spot
size applied to tattoos over the treatment period.

  * •

Greater average fluence applied to a tattoo is related to increased likelihood
of complication occurrence.

  * •

Greater average laser wavelength over a treatment period is both
parametrically and non-parametrically related to increased likelihood of
complication occurrence. Shorter wavelengths may be emulating intense
radiation causing complications associated with skin discoloration.

  * •

According to our logistic regression model, average wavelength, standard
deviation in laser spot size, and overall change of laser frequency applied to
a tattoo through a given treatment sequence are related to increased
likelihood of complication occurrence.

It is also worth noting the response-factor relationships that were found to
not be statistically significant with respect to the current literature
(again, with brief explanations and implications):

  * •

The average number of days between appointments is not statistically related
to complication occurrence. The number of days between appointments in the
data is sufficiently long for skin recovery before a follow-up appointment.

  * •

Tattoo age and patient age is not statistically related to the variation in
complication rates. This possibly suggests that current clinician approaches
to tattoo removal take into account tattoo and patient age, since the
application of certain settings used by clinicians vary between different ages
at a statistically significant level.

  * •

Between groups of tattoos that experienced complications and those that did
not, the variation in average laser fluence, overall change in laser fluence
between the first and last appointment, average spot size, and overall change
in laser spot size between the first and last appointment were not found to be
distinct at a statistically significant level.

  * •

Patient sex, ethnicity, Fitzpatrick score, and being above the median age are
not statistically related to complication occurrence within a recorded
treatment sequence, independent of whether a patient completed the entire
treatment process.

Finally, note the discrepancy of our results between the two machine learning
models. Although our boosted decision tree procedure lowly ranks a
statistically significant variable—according to our logistic regression
model—such as “Mean Difference Laser Frequency,” this result is not an
inconsistent finding. The two models essentially offer distinct
interpretations of the data.

The logistic regression model measures that in a total interaction of all ten
variables, average wavelength, standard deviation in spot size, and average
variation of laser frequency between the first appointment and pre-arrival of
a recorded tattoo complication are statistically significant factors for the
tattoo removal practitioner to be aware of in the treatment process. On the
other hand, the aggregated boosted decision tree model ranks variables by
observing how our predictions of complication occurrence/non-occurrence are
affected if they are not considered in the interaction. For example, excluding
average wavelength in the interaction of treatment parameters will most
frequently negatively affect classification of complication
occurrence/nonoccurrence in our data compared to the exclusion of average
frequency which will be of least consequence in our predictions.

Limitations worth noting with respect to our findings are that these results
are specific to our given dataset, which has a fair share of missing values
without any data imputation having been performed, and we do avoid making any
full clinical interpretations of our results as well. Yet, our data is
relatively abundant and detailed, and we’ve applied many interesting
techniques along with applying this ad hoc aggregated variable ranking
procedure using gradient boosted decision trees, which may prove to be a
theoretically well-justified and statistically powerful way of producing ranks
of variable importance.

### Acknowledgements

JH, DN and RY are grateful to Jessica Bogner at Homeboy Industries, Los
Angeles, and Dr. Jo Marie Reilly at Keck Medicine of USC, for providing data
and guidance with this project. The authors were partially supported by NSF
BIGDATA #1740325 and NSF DMS #2011140.

## References

  * [1] W Bäumler and KT Weiß.  Laser assisted tattoo removal–state of the art and new developments.  Photochemical & Photobiological Sciences, 18(2):349–358, 2019\.
  * [2] David Roxbee Cox.  Principles of statistical inference.  Cambridge university press, 2006.
  * [3] Scott A Czepiel.  Maximum likelihood estimation of logistic regression models: theory and implementation.  Available at czep. net/stat/mlelr. pdf, pages 1825252548–1564645290, 2002.
  * [4] C Montgomery Douglas.  Design and analysis of experiments.  John Wiley and Sons Inc, 2001.
  * [5] Bradley Efron and Robert J Tibshirani.  An introduction to the bootstrap.  CRC press, 1994.
  * [6] Thomas B Fitzpatrick.  Soleil et peau.  J Med Esthet, 2:33–34, 1975.
  * [7] Phillip Good.  Permutation tests: a practical guide to resampling methods for testing hypotheses.  Springer Science & Business Media, 2013.
  * [8] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  The elements of statistical learning: data mining, inference, and prediction.  Springer Science & Business Media, 2009.
  * [9] Stephanie GY Ho and Chee Leok Goh.  Laser tattoo removal: a clinical update.  Journal of cutaneous and aesthetic surgery, 8(1):9, 2015.
  * [10] Gareth James, Daniela Witten, Trevor Hastie, and Robert Tibshirani.  An introduction to statistical learning, volume 112.  Springer, 2013.
  * [11] Clyde Young Kramer.  Extension of multiple range tests to group means with unequal numbers of replications.  Biometrics, 12(3):307–310, 1956.
  * [12] William H Kruskal and W Allen Wallis.  Use of ranks in one-criterion variance analysis.  Journal of the American statistical Association, 47(260):583–621, 1952.
  * [13] Ivan Kurniadi, Farida Tabri, Asnawi Madjid, Anis Irawan Anwar, and Widya Widita.  Laser tattoo removal: Fundamental principles and practical approach.  Dermatologic Therapy, page e14418, 2020.
  * [14] Henry B Mann and Donald R Whitney.  On a test of whether one of two random variables is stochastically larger than the other.  The annals of mathematical statistics, pages 50–60, 1947.
  * [15] Mary L McHugh.  The chi-square test of independence.  Biochemia medica, 23(2):143–149, 2013.
  * [16] Jørgen Serup and Wolfgang Bäumler.  Guide to treatment of tattoo complications and tattoo removal.  In Diagnosis and Therapy of Tattoo Complications, volume 52, pages 132–138. Karger Publishers, 2017.
  * [17] Robert H Shumway, David S Stoffer, and David S Stoffer.  Time series analysis and its applications, volume 3.  Springer, 2000.
  * [18] R Core Team et al.  R: A language and environment for statistical computing.  2013\.
  * [19] John W Tukey.  Comparing individual means in the analysis of variance.  Biometrics, pages 99–114, 1949.
  * [20] Ronald E Walpole, Raymond H Myers, Sharon L Myers, and Keying Ye.  Probability and statistics for engineers and scientists, volume 5.  Macmillan New York, 1993.
  * [21] Bernard Lewis Welch.  On the comparison of several mean values: an alternative approach.  Biometrika, 38(3/4):330–336, 1951.


###### Abstract

Currently, space debris represents a threat for satellites and space-based
operations, both in-orbit and during the launching process. The yearly
increase in space debris represents a serious concern to major space agencies
leading to the development of dedicated space programs to deal with this
issue. Ground-based radars can detect Earth orbiting debris down to a few
square centimeters and therefore constitute a major building block of a space
debris monitoring system. New radar sensors are required in Europe to enhance
capabilities and availability of its small radar network capable of tracking
and surveying space objects and to respond to the debris increase expected
from the New Space economy activities. This article presents ATLAS, a new
tracking radar system for debris detection located in Portugal. It starts by
an extensive technical description of all the system components followed by a
study that estimates its future performance. A section dedicated to waveform
design is also presented, since the system allows the usage of several types
of pulse modulation schemes such as LFM and phase coded modulations while
enabling the development and testing of more advanced ones. By presenting an
architecture that is highly modular with fully digital signal processing,
ATLAS establishes a platform for fast and easy development, research and
innovation. The system follows the use of Commercial-Off-The-Shelf
technologies and Open Systems which is unique among current radar systems.

###### keywords:

radar, sst, ssa, space debris, debris, leo, tracking, doppler, eusst, ATLAS

xx 1 5 Received: date; Accepted: date; Published: date Development of the
first Portuguese radar tracking sensor for Space Debris João Pandeirada 1,2,
Miguel Bergano 1,3, João Neves 4, Paulo Marques 5, Domingos Barbosa 1, Bruno
Coelho 1 and Valério Ribeiro 1,2 João Pandeirada, Miguel Bergano and João
Neves

##  1 Introduction

The use of sensors (which include telescopes, lasers and radars) has been a
growing area of interest for monitoring the space environment around Earth.
Due to the growing number of space debris, there is a need to predict and
adjust the orbits of the more than 2000 active satellites Statista , avoiding
collisions with other inactive satellites or debris, to guarantee their long-
term operation and investment.

These actions are core activities of the Space Surveillance and Tracking (SST)
concept or, in a broader sense, are framed within the Space Situational
Awareness (SSA) domain, notwithstanding the different concepts acknowledged
between a USA concept, strongly related to the Defense, or the European
concept, strongly related to a dual-use or purely civil use compliant with the
European Union (EU) or the European Space Agency (ESA) guidelines.

The Earth near-space environment has been getting polluted with space debris
since the onset of the space exploration era in 1957. The envisaged space
launching actions to put more satellites in orbit including the announced mega
constellations or broadband internet satellites, e.g. Starlink Space X and One
Web, is a growing concern which drives the need to create, maintain and
improve accurate sensor networks to collectively contribute to a safe use of
space. The safety and security of space assets is an important matter of
national interests and central to international treaties and regulations,
namely those agreed to in the framework of the United Nations, at the
Committee on the Peaceful Uses of Outer Space United Nations Office for Outer
Space Affairs .

Sensors networks comprise both ground based and space based assets, although
only a few nations can afford to have space based ones. Regarding the ground
based sensors, these are usually divided in three main groups: telescopes,
radars and lasers.

Radar usage for SST and SSA purposes is increasing due to its capacity for
monitoring the space environment during day and night, and insensitivity both
to light pollution and to weather conditions. A tracking radar distinguishes
itself not only by recognizing the presence of a target and determining a
target location in range and in angle coordinates, but also because it is able
to follow the target and observe it over some time range thus improving the
information accuracy about its trajectory Skolnik (2001).

The principle of tracking radar is to use the angular error signal, defined by
the difference between the target direction and the reference direction given
by the axis of the antenna, to adjust the antenna’s pointing direction. The
tracking radar attempts to position the antenna with zero angular error (i.e.,
to locate the target along the reference direction) Chen (2004) .

The first developments of tracking radars date back from World War II Galati
(2016) but their biggest developments date from the Cold-War Ingwersen and
Lemnios (2000). During those times, both the USA and the Soviet Union built
large networks of ground-based tracking facilities that collected data on man-
made objects in Earth orbit. They consisted, primarily, of large tracking
radars and optical telescopes, often dual-purposed with other military
missions, such as ballistic missile warning and tracking. The main focus of
space surveillance at that time was on determining the precise orbit of man-
made objects around the Earth Schrogl et al. (2015). Later, they were improved
to make them useful for Combat Weapon Systems, such as the Thales STIR-
Tracking and Illumination radar Thales Group .

The tracking of distant objects in the higher layers of the atmosphere, in the
lower orbits around Earth or even far distant objects throughout the galaxy,
continued to develop up to complex phased array systems such as the Square
Kilometer Array (SKA) radio telescope SKA  or the brand new US Space Fence,
costing several million dollars eoPortal , to name a few. Tracking radars are
now suitable for Low-Earth Orbit (LEO) tracking or space object catalogue
accuracy improvement due to electronic components developments and increased
digital processing capability.

Referring to existing operational tracking radars for in-orbit object
detection (and disregarding the ones based on phased arrays, which usually
also have surveillance capability), these can be found all over the world and
with different capabilities.

On the US side, the development of these radars has a long history with the
development of several radars, most of them for the Ballistic Missile Defense
program, supported by DARPA (Defense Advanced Research Projects Agency), which
were closely related to the detection of objects in orbit Ingwersen and
Lemnios (2000). Notwithstanding the huge variety of radars that were developed
and used (some of them have been decommissioned), two of the most interesting
and relevant radars used for Space Debris detection, tracking and
characterization are the Haystack and the HAX (Haystack Auxiliary Radar)
radars. The HAX radar has a 12.2m antenna, operating at 16.7GHz and is much
smaller than Haystack which has a 36.6m antenna operating at 10GHz Stokely et
al. (2006) , both in a Cassegrain configuration. HAX had twice the bandwidth
of Haystack, being capable of producing correspondingly sharper images of
near-earth orbit objects. The cost of building HAX, which became operational
in 1994, was reduced by using a refurbished antenna and by sharing signal
processing and data processing systems with the Haystack Long Range Imaging
Radar (LRIR) (means that the two radars could not operate simultaneously). The
Haystack radar has been operational since 1964, but from 2010 to 2014 it
underwent several upgrades becoming the Haystack Ultra-Wideband Satellite
Imaging Radar (HUSIR), able to also operate in the W-band with a bandwidth up
to 8GHz, which allowed the HAX to detect objects down to 3cm and the HUSIR
down to 3mm at an altitude up to 1000km Murray et al. (2019).

The Space Debris networks around the world includes facilities from other
major players, such as China Du (2017) and Russia and other spacefaring
nations such as India Prasad (2005) Adimurthy and Ganeshan (2006), Iran,
Brazil Froehlich et al. (2020) Nogueira et al. (2012), Japan or North Korea.
However, most of those nations projects’ and infrastructures are linked to the
military and involved in deep secrecy. It’s highly probable that radio-
telescopes and other related ground-based sensors support missions including
intelligence collection, counterspace targeting, ballistic missile early
warning, spaceflight safety, satellite anomaly resolution, and also space
debris tracking DEFENSE INTELLIGENCE AGENCY (2019) (DIA).

In the European continent the development of tracking radars for orbit objects
detection has evolved for some time and nowadays there are around 50 different
radio telescopes. Just to mention some which have played a relevant role in
space debris and satellite detection, it is worth mentioning the German
Tracking and Imaging Radar (TIRA) system, with a 34 meter antenna, located
near Bonn. In a monostatic configuration it is able to detect objects of
diameters down to 2 cm at 1000 km, but in a bistatic configuration, in synergy
with the Effelsberg Radio Telescope of the Max Planck Institute for Radio
Astronomy, which has a 100 meter dish antenna, the system can detect objects
down to 1 cm diameter Mehrholz (1997) Mehrholz et al. (1999). The TIRA radar
has been one of the most useful radars in Europe in the SST context and its
accuracy and imaging capabilities, along its continued upgrades, made it a
valuable asset. The Sardinia Radio Telescope (SRT) is an Italian radio
telescope facility with recent radar capabilities with a 64 meter antenna
located in the Sardinia island, operating only since 2012. Compared to
equivalently sized radars, it has a fairly good speed in azimuth and elevation
tracking (0.85º/s and 0.5º/s, respectively) and is able to operate from 300MHz
up to 15 GHz. It is intended to work in a bistatic configuration, using the
Flight Terminator System (FTS), about 40 km apart from the first, which
configuration is named BIRALET (BIstatic RAdar for LEo Tracking) Muntoni et
al. (2017). BIRALET has been used in the context of the European SST for LEO
detection, but its capabilities, namely the SRT antenna, envisage the
possibility to go beyond this orbit.

On the French side, one of the countries in Europe that has more investments
in the space segment, there have been several developments in radars, but the
most relevant and interesting of them are in the military side. The French
space surveillance network comprises a set of three tracking radars named
SATAM, in three different locations (one of it with deployable capability).
This configuration is not primarily dedicated to space surveillance, but is
used for space events, detection of risk collision and atmospheric reentries,
being its speed (about 40º/s) an incomparable added-value, notwithstanding the
lower dimension of the dishes, compared with the previous ones.

The European Space Surveillance and Tracking (EUSST) is a support framework
setup by the European Commission to create an autonomous monitoring network
supporting the detection and tracking of space debris and issue alerts when
evasive actions may be necessary. EUSST aims to establish an SST capability at
the European level with an appropriate level of European autonomy European
Union (2014),EU SST . Portugal is a member of the EUSST program and is
developing capabilities both in optical and radar sensors.

This paper describes the development of the first Portuguese tracking sensor
for space debris and is organized as follows. Section 1 explains the necessity
of monitoring the space environment and how nations worldwide have tackled
this problem throughout the years. It also provides examples of already
deployed tracking radars with similar functions to the current one in
development. Section 2 introduces ATLAS (rAdio TeLescope pAmpilhosa Serra) and
provides a detailed technical description of its components: antenna,
transmitter, receiver, clock generation, data acquisition and controllers.
Section 3 consists of simulation studies that use the technical
characteristics of the system in order to estimate some performance metrics:
minimum detectable target size, minimum elevation, number of expected
observable objects and maximum simultaneous trackable targets. Section 4
explains the importance of waveform design on radar systems and provides some
examples of waveforms that can be used by ATLAS. Finally, Section 5 concludes
the article and reveals the next steps on the realization of this project.

##  2 Implemented System

Figure 1: ATLAS system block diagram

In order to support and consolidate the pooling of national resources to
prepare for the establishment of a SST sensor network, Instituto de
Telecomunicações initiated in 2019 a major upgrade of its Cassegrain Antenna
located in Pampilhosa da Serra, Portugal. Attached to the antenna will be a
complete ground-based radar system, named ATLAS, operating in the C-band (5.56
GHz).

ATLAS is a monostatic pulse radar using solid state power amplifiers (GaN)
with a peak output power of 5 kW for tracking debris objects in LEO. Regarding
the receiver side, the system is fully coherent with detection and processing
fully in the digital domain with a bandwidth of 50 MHz and with the capacity
to detect Doppler velocities up to 10.79 km/s. The architecture, depicted in
Figure 1, is highly inspired in the HAX/Haystack system Eshbaugh et al. (2014)
developed by NASA , which is a very successful ground-based radar used for
debris detection and tracking. The system follows the use of Commercial-Off-
The-Shelf (COTS) technologies and Open Systems (OS) which enables a major
downstream cost reduction in the development and maintenance of the system
Rejto (2000). Table 1 provides a comparison between ATLAS and other tracking
radars currently used in SSA activities.

Table 1: Comparison between ATLAS and other tracking radars | ATLAS | TIRA | BIRALET | HAX | Haystack
---|---|---|---|---|---
Operating Frequency | 5.56 GHz | 1.333 GHz | 410 MHz | 16.7 GHz | 10 GHz
Peak power | 5 kW | 1 MW | 4 kW | 50 kW | 250 kW
Waveform type | Pulsed | Pulsed | Continuous | Pulsed | Pulsed
Antenna Gain | 46 dB | 49.7 dB | 13 dB (TX), 47 dB (RX) | 63.64 dB | 67.23 dB
Antenna Beamwidth | 0.73º | 0.50º | 30º (TX), 0.8º (RX) | 0.1º | 0.058º
Receiver Bandwidth | 80 MHz | 250 kHz | 5 MHz | 1 MHz | 1 MHz
Topology | Monostatic | Monostatic | Bistatic | Monostatic | Monostatic

###  2.1 Radar site and Antenna

We obtained a Vertex RSI C band, 9 metre Cassegrain Antenna originally used
for Satellite Communications in the Azores Terceira Island, Portugal. The
antenna was moved to Pampilhosa da Serra, Portugal, a site with a low Radio
Frequency Interference (RFI) environment Fonseca et al. (2006). New
foundations, electrical power and lightning protection were installed at this
new site. The mounting pedestal underwent a major adaptation in order to
comply with a fast, continuous azimuth antenna rotation. Originally intended
for radio astronomy surveys, this facility represents a step forward in our
technological development and evolution of high-performance space radar
instrumentation systems.

The major improvements have given the station space telemetry and tracking
capabilities. These improvements will represent one of Portugal’s main support
pilot facilities for space monitoring and will provide support to current and
future space missions. The upgrade approach manages to automatically, and
quickly, attain the best antenna steering strategy to successfully track the
objects.

The data output capabilities have been upgraded for wide bandwidth. These
capabilities are now available to users on site and remotely via a high
throughput connection. Extensive user support will be provided for these new
facilities, including assistance, monitoring and control, antenna interfaces
and output interfaces. Operation and data reduction is also available for
object tracking.

The antenna is depicted in Figure 2 and specifications are detailed in Table
2.

Table 2: Antenna specifications

Optical Configuration | Cassegrain
---|---
Mount configuration | Alt-Azimuthal
Primary aperture | 9.0 m
Primary depth | 1.53
Secondary aperture | 1.17
F/D | 0.368
Surface rms (static) | 0.5 mm
Azimuth Travel | 0º to 360º
Azimuth Travel Rate | 4 rpm
Elevation Travel | 35º to 90º
Elevation Travel Rate | 2º/s
Maximum Operational Wind (5GHz) | 25 km/h
Survival Winds | 150 km/h
Reflector Weight | 1900 kg
Pedestal Weight | 3500 Kg
Foundation Size | 4.5 x 4.5 x 0.75 m
Concrete Volume | 12 ms
Beamwidth @ 5 GHz | 44 arcmin
Beamwidth @ 10 GHz | 22 arcmin
Pointing Accuracy (wind limited) | <1/20 beamwidth
Gain @ 5 GHz ($G$) | 46 dB
G/T (5 GHz) | 39 dB/K

Figure 2: ATLAS Antenna

###  2.2 Transmitter

Active radars require Power Amplifiers (PAs), which should be small, efficient
and low-cost. We implemented a C-band PA centered at 5.56GHz with an output
peak power of 5 kW (67dBm) from a 30dBm drive with a Power Added Efficiency
(PAE) of 67 %. GaN technology is suitable for these requirements and can do so
with high efficiency. The PA was implemented by readily available COTS
technology. Compact Radio Frequency (RF) power sources now exist with the
recent release of high-power solid-state RF amplifiers such as the
Wolfspeed/Cree CGHV59350 high-electron-mobility transistors (HEMTs) Cree ()
(2018). These HEMTs are capable of 500 W of RF power each. Our design follows
the traditional approach of combining a large number of solid-state RF
amplifiers into a high-power all-solid-state RF system Simons et al. (2019).
The transmitter, depicted on Figure 1 in blue, is responsible for outputting
the RF signal. All the relevant operating parameters are summarized in Table
3. The signal generation starts with a frequency synthesizer at 5.160 GHz
which is fed to a frequency translator. Along with the 5.160 GHz signal, the
frequency translator receives an intermediate frequency of 400 MHz from the
reference module and generates the desired 5.560 GHz signal. This signal is
fed into a PIN modulator (PIN diode based) responsible for modulating the
5.560 GHz carrier in amplitude. The PIN modulator is also connected to the
controller board in order to receive the user-defined pulse shape. The
modulated signal goes to a driver, to achieve a 30dBm minimum signal to be
provided to the PA. We selected the previously mentioned Cree HEMTs at 5–6 GHz
because these HEMTs provide the highest available RF power as required
(67dBm). The target performance requirements for the amplifier are listed in
Table 3.

Table 3: Transmitter Operating Parameters Peak power | 5 kW
---|---
Transmitter Frequency ($1/\lambda$) | 5.56 GHz
Waveform | Arbitrary amplitude modulation
Max. Pulse length ($\tau$) | 10 s
Phase Noise | -91.3dBc[Hz] @ 100 kHz
Modulator | Modified D195 General Microwave
PA Transistors | CGHV59070 Wolfspeed Cree , CGHV59350 Cree () (2018)

###  2.3 Receiver

The antenna is equipped with a corrugated feed horn, followed by a polarizer
and by an orthomode transducer (OMT) to carry the left and right circular
polarizations (LHCP for Tx and RHCP for Rx). A waveguide to coaxial transition
performs the connection to the receiver. The receiver has a limiter and switch
right at the input that handles an excess peak power of 20 W for about 10 ms
and blocks the signal during transmission with an attenuation of 30 dB.

The receiver, depicted on Figure 1 in purple, includes an analog chain that is
responsible for amplifying the signals (target echoes) to the required digital
entry levels of the acquisition board (yellow). The receiver starts with an
low noise amplifier (LNA), internally developed, using GaAs technology,
presenting a noise factor of 0,7 dB at 5.56GHz. Next, there is a local
oscillator at 5.160 GHz to convert it down to a 400 MHz intermediate frequency
(IF). The local oscillator is provided by the same reference signal from the
transmitter, guaranteeing a coherent conversion. The IF then passes through an
IF filter and a digital variable gain amplifier (DVGA) which is connected to
the controller board.

Finally, the I/Q detector converts the signal to base-band (or 400 kHz) and
outputs the I and Q components of the echo signal. The conversion to 400 kHz
is used for accommodating Doppler shifts, providing velocity measurements up
to 10.79 km/s. Table 4 summarizes all relevant parameters of the receiver
chain.

Table 4: Receiver Operating Parameters IF | 400 MHz
---|---
LNA Noise figure (<15ºC) ($N_{f}$) | 0.7 dB
Receiver Temperature ($T$) | 15ºC
IF Filter BW ($B$) | 80 MHz @ -3dB
IQ Detector output BW | 50 MHz

###  2.4 Reference

The reference module, depicted on Figure 1 in orange, generates reference
signals at 100, 400 and 400.4 MHz. The main reference oscillator is a 5 MHz
oven-controlled crystal oscillator (OCXO) from MTI MTI-Milliren Technologies,
Inc. . The 5 MHz reference locks a pair of 100 MHz XPLO’s that are then
multiplied by 4 to generate the 400 and 400.4 MHz. This module is connected to
the controller board in order to switch between IQ detection at base band or
with 400 KHz offset.

###  2.5 Controller

The controller, depicted on Figure 1 in green, is responsible for the digital
control of the whole system. It is composed of a controller board, a pulse
generator and an arbitrary waveform generator (AWG).

The board establishes a telnet connection with the host PC via Ethernet (ETH)
where the user can program the radar using a command line interface (CLI). The
CLI can be used to define a great range of radar parameters such as the pulse
length, shape, time between pulses, number of pulses and a lot of other
internal parameters such as gains and delays. It is also possible to display
several system status parameters such as pulse configurations and
temperatures. Table 5 shows all the variables that can be monitored and
configured.

The pulse generator is responsible for timing all the components of the
system. It enables/disables the transmitting and receiving chains and triggers
the arbitrary wave generator as well as the acquisition board. This enables
the system to be fully coherent and synchronized from signal generation to
acquisition.

Finally, the AWG enables the design of arbitrary amplitude modulated waveforms
in the digital domain that can be uploaded to the system.

Table 5: Controller board parameters Waveform resolution | 10 ns
---|---
Pulse repetition frequency ($f_{p}$) | 10 MHz (max)
Duty Cycle | 10% (max)
Number of pulses for integration | Variable
Temperature Monitoring | Tx, Rx, LNC, Driver, PA, Coupler
System Monitoring |  | Pulse shape file name, Tx/Rx gain,
---
No. Pulses, PRF, Pulse length …

###  2.6 Data acquisition

The data acquisition, depicted on Figure 1 in yellow, contains the acquisition
board (AQS) and a single board computer (SBC). The acquisition board is
triggered by the controller board and is responsible for converting the I/Q
signal components to the digital domain, by employing two 16 bit ADC’s at a
maximum sample rate of 100 MS/s. The computer (MIO 2360 N-32A1E by Advantech
Advantech ) receives the I/Q signal in the digital domain and performs the
necessary digital signal processing. Communication between the AQS and the SBC
is done via full speed USB 3.0 (100 MS/s) and a Fast/Slow mode buffer is
selectable, where a compromise between sample resolution and sample rate can
be chosen.

##  3 System expected capabilities

Given the technical specifications of the current system, it is possible to
make estimates on how it will perform. These estimates are not only useful for
predicting the systems applicability but also to identify its limitations and
guide to future upgrades.

###  3.1 Minimum detectable target size

The radar range equation represents the physical dependencies of the
transmitted power and is used to obtain the power in the receiving antenna
Skolnik (2008). Using the radar equation, we can obtain the dependence of the
signal to noise ratio (SNR) with the radar specifications:

$\text{SNR}=\frac{P_{av}G^{2}\lambda^{2}\sigma ne(n)F^{4}}{4\pi^{3}\tau
f_{p}R^{4}N_{F}kTBL_{s}}$ (1)

where:

  * •

$P_{av}$ is the average transmitted power ($W$)

  * •

$G$ is the antenna gain

  * •

$\lambda$ is the operating wavelength ($m$)

  * •

$\sigma$ is the target radar cross section (RCS) ($m^{2}$)

  * •

$n$ is the number of integrated pulses

  * •

$e(n)$ is the integration efficiency

  * •

$F$ accounts for all the propagation effects

  * •

$\tau$ is the pulse width ($s$)

  * •

$f_{p}$ is the pulse repetition frequency ($Hz$)

  * •

$R$ is the distance to the object ($m$)

  * •

$N_{F}$ is the noise factor of the receiver

  * •

$k$ is the Boltzmann constant ($\frac{J}{K}$)

  * •

$T$ is the receiver temperature ($K$)

  * •

$B$ is the bandwidth of the receiver ($Hz$)

  * •

$L_{s}$ accounts for all system losses

Some of the radar specifications are fixed and have already been shown during
this article. Regarding propagation losses it is known that for frequencies
below 10 GHz losses due to atmospheric absorption may be neglected
International Telecommunication Union , for that reason we neglected for now
the $F$ factor. In future work we will include a detailed modelling of the
atmospheric propagation effects including tropospheric absorption. Accounting
for all the losses on the waveguide transition from the antenna to the
polarizer/orthomode transducer (OMT)/switch and the coaxial transition from
the transmitter/receiver to the switch, we obtained $L_{s}=0.6$ dB.

By changing the pulse length, pulse repetition frequency and the number of
integrated pulses, it is possible to establish a desired SNR for an
established target RCS. By doing a statistical analysis on the mean elevation
speed of the debris obtained from the Space-Track catalogue (explained in 3.3)
and considering our antenna beamwidth, we can safely define 20 pulses for
integration. Since the radar is fully coherent with integration in the digital
domain, we can assume a coherent integration with full efficiency ($e(n)=1$).
In order to obtain the desired SNR, we defined $\tau=3.3$ ms and $f_{p}=30$
Hz, obtaining the following performance results:

  * •

Maximum Unambiguous Range: 5000 km

  * •

Maximum Unambiguous Velocity: 0.4 m/s

  * •

SNR for a 1 m2 RCS at 103 km: 39.55 dB

  * •

SNR for a 10 cm2 RCS at 103 km: 9.55 dB

We can conclude that we can measure targets up to 1000 km in distance with a
minimum RCS of 10 cm2.

Since debris objects in LEO have velocities of around 7.8 km/s, obtaining a
maximum unambiguous velocity of 0.4 m/s results in inability to measure the
Doppler shift. This is due to the fact that we need a low pulse repetition
frequency (PRF) for maintaining range unambiguity. It is still possible to
measure velocities without Doppler shift by taking several range measurements
and inferring velocity from the measurements, which in fact is what is done in
an SST system, by processing several individual tracks of the same object.

In order to take advantage of the long idle times between the transmitted and
received pulses in a low PRF system, we can use a process called pulse
interleaving. Pulse Interleaving is a technique that inserts carefully crafted
pulses in the idle time of other pulses. One example is the emission of
different waveforms with low cross correlation between them, such as
Orthogonal Frequency Division Multiplexing (OFDM) signals, in order to
suppress range ambiguity Riché et al. (2014).

###  3.2 Elevation for initial and final target acquisition

The location of the radar and its surroundings need to be taken into account
for object tracking. ATLAS is located at (40.18ºN,7.87ºW) and is surrounded by
mountains, which makes the horizon altitude profile vary in some directions.
In these directions, this decreases the arc length of a debris passage thereby
reducing its potential tracking time. We thus defined that the minimum
elevation for initial and final target acquisition as 30º.

###  3.3 Number of expected observable debris objects

ATLAS was designed to track space objects in LEO orbits below 1000 km of
altitude, however, several objects in this range might not be observable due
to several constraints:

  * •

Objects with an RCS below the minimum detectable threshold;

  * •

Orbits with an elevation range below the minimum elevation required for
detection and data acquisition;

  * •

Orbital speeds exceeding the maximum antenna tracking speed;

In order to make an estimate on the number of objects trackable by ATLAS, we
retrieved the latest two-line element (TLE) files (which contain the orbital
elements) from all the objects categorized as debris with an apogee less than
1000km from the Space-Track public catalogue USSPACECOM . With the orbital
elements, we predicted the orbit of each object for the next 7 days (counting
from the epoch of the TLE file) with an SGP4 based algorithm Satellogic SA .
Then, we intercepted the predicted orbits with the visible portion of the sky
for ATLAS, which is obtained by the information provided in Section 3.2.
Objects that do not intercept the visible sky are discarded.

Since the catalogue does not provide detailed information about the objects
RCS, it is not possible to filter out by RCS. All the objects were considered
regardless of RCS and the filtering was done by minimum elevation and maximum
velocity. Figure 3 depicts the number of expected observable objects during a
complete day.

Figure 3: Number of expected observable debris objects during a 24 hour
timespan (regardless of RCS).

###  3.4 Maximum simultaneous number of trackable targets

Since the antenna has a narrow beamwidth (see Table 2), the amount of sky
visible at a specific moment in time is limited and provides an idea on the
number of objects illuminated by the radar simultaneously. Figure 4 shows the
maximum number of simultaneous trackable targets over the course of the day
for a beampark configuration. Figure 3 and 4 share the same shape of the
distribution because, if we assume that the number of objects at each bin in
Figure 3 is uniformly distributed over the visible sky, the number of
simultaneous objects is a fraction of that number, thus the shape is retained.

Figure 4: Maximum simultaneous number of trackable targets during a 24 hour
timespan.

##  4 Waveform design

Designing a proper waveform for the radar pulse is of utmost importance. By
manipulating the transmitted pulses in terms of amplitude, frequency and phase
it is possible to create waveforms that maximize the range and velocity
resolution of the echo signal. In order to evaluate the characteristics of a
waveform one can use the ambiguity function (AF). The AF represents the time
response of the matched filter when the signal received is affected by a delay
$d$ and a Doppler shift $v$ relative to the values expected by the filter and
is given by:

$|X(d,v)|=\bigg{|}\int_{-\infty}^{\infty}u(t)u^{*}(t+d)e^{j2\pi vt}dt\bigg{|}$
(2)

where $u$ is the complex envelope of the signal N.Levanon and E.Mozeson
(2004a). By analyzing the zero delay and zero Doppler cuts on the AF it is
possible to define the range and velocity resolution. Since ATLAS will work
with a low pulse repetition frequency when detecting debris, waveform design
focuses on maximizing range resolution instead of velocity.

Matched filtering is the most basic form of signal processing commonly used in
radar. It consists in the convolution of the received signal with a filter
that is matched with the emitted signal. Matching will result in the maximum
attainable SNR at the output of the filter when the signal to which it was
matched, after addition of white noise, is passed through it N.Levanon and
E.Mozeson (2004b).

As defined in Section 3.1, ATLAS will use pulses of 3.3 ms with a pulse
repetition frequency of 30 Hz in order to maintain the desired SNR. The
waveform used can be arbitrary, within the hardware limits described
previously and as long as those power requirements are met.

###  4.1 Rectangular Pulse

Figure 5 shows a rectangular pulse of 3.3 ms repeated at 30 Hz and the zero
Doppler cut. Using a rectangular pulse results in a range resolution of 500 km
which is unacceptable when detecting objects at a maximum of 1000 km.

Figure 5: Rectangular pulse (left) and zero cut of AF (right).

ATLAS enables the design of arbitrary amplitude modulated waveforms in the
digital domain with a resolution of 10 ns. This enables the design of
waveforms with high autocorrelation properties that can be used to increase
resolution. This signal processing technique is usually called pulse
compression Skolnik (2008).

###  4.2 Phase-Coded Pulse (Barker Code)

The Barker code is a famous binary code used in pulse compression and consists
in binary sequences that guarantee that the peak-to-peak sidelobe ratio of the
autocorrelation is $M$ where $M$ is the size of the Barker code. Another
important feature of the Barker code is that the range resolution increases by
$M$ in relation to the rectangular pulse. We can further increase this
resolution by using nested Barker codes N.Levanon and E.Mozeson (2004).

Figure 6 shows a Nested Barker code composed of two 13 element codes and the
corresponding zero Doppler cut. For visualization purposes, only a small
fraction of the code inside of the 3.3 ms pulse is shown. As we can see, the
pulse is compressed by a factor of 169, resulting in a range resolution of
approximately 2.96 km which is a massive improvement from the rectangular
pulse.

Figure 6: Part of the Nested Barker Code (left) and zero cut of AF (right).

Longer binary codes with good auto correlation properties are available which
enable even higher compression factors Dzvonkovskaya and Rohling (2008).

###  4.3 Linear Frequency-Modulated Pulse

Linear frequency modulation (LFM), also known as chirp, is probably the most
popular waveform used for pulse compression. It consists in sweeping the
carrier wave by the frequency band $B$ during the pulse duration $T$. The
complex envelope of a linear chirp is given by N.Levanon and E.Mozeson (2004):

$u(t)=\frac{1}{\sqrt{(}T)}rect(\frac{t}{T})e^{j\pi kt^{2}}$ (3)

where $k=\pm\frac{B}{T}$ for a linear chirp and the signal indicates an
increasing/decreasing sweep. The time-bandwidth product of the signal $B\times
T$ is the compression factor in relation to the rectangular pulse, that is, by
increasing $B$ the range resolution improves. Since ATLAS currently only
supports amplitude modulation of the carrier, it is not possible to generate
traditional chirp pulses because it requires phase modulation. In order to
overcome this limitation, we developed an AM Chirp.

The AM Chirp consists in designing a chirp signal at a much lower frequency
and using it for modulating the carrier in amplitude. The amplitude modulation
is given by:

$y(t)=A(t)\cos(2\pi f_{c}t)$ (4)

where $A(t)$ is the modulating signal and $f_{c}$ is the carrier frequency. In
the case of the AM chirp, $A(t)$ is given by:

$A(t)=\alpha+\beta\cos[2\pi(ct+f_{0})t],~{}c=\frac{f_{1}-f_{0}}{T}$ (5)

where $\alpha$ and $\beta$ define the modulation parameters, $f_{0}$ and
$f_{1}$ are the starting and final frequencies and $T$ is the pulse duration.

Figure 7 illustrates a simple case of AM chirp. The left side depicts a 1 kHz
carrier modulated in amplitude by a chirp signal with $\alpha=0.7$,
$\beta=0.3$, $f_{1}=0$ and $f_{1}=50$ Hz. After carrier demodulation,
filtering and signal reconstruction in the digital domain, the chirp waveform
is successfully recovered as depicted on the right side of Figure 7.

Figure 7: AM Chirp (left) and reconstructed chirp (right).

As said in Section 2, ATLAS has a waveform resolution of 10 ns and a receiver
bandwidth of 80 MHz, which enables generation of AM chirps with high
compression factors. As an example, Figure 8 shows a linear chirp with a
compression factor of 1000 inside of the 3.3 ms rectangular pulse. For
visualization purposes, only a small fraction of the waveform inside of the
pulse is shown.

A compression factor of 1000 was used to maintain a reasonable simulation time
while showing the benefits of using this waveform for pulse compression. As
can be seen on the right side of Figure 8, the range resolution now is
approximately 500 m which corresponds to an error of 0.05% at 1000 km of
altitude.

Figure 8: Part of the chirp waveform (left) and zero cut of AF (right).

This waveform requires a bandwidth of 300 kHz which is only a small fraction
of the capabilities of the system to generate and receive it. This indicates
it is possible to synthesize waveforms with even higher compression factors.

###  4.4 Advanced Waveform Design

In the particular scenario where certain selected man-made targets are of
interest, the transmitted waveform can be further optimized Marques (2017).
The basic idea is to change the signal used as reference for detection in such
a way that it is adapted to the predominant scatterers of the targets,
resulting in a better matched filtering which ultimately leads to better SNR
and pulse compression.

ATLAS also permits the testing of waveforms typically used on Noise Radar
Technology K. Savci et al (2020), being able to transmit a virtually infinite
set of waveforms produced via a set of random realizations, with the real-time
configuration of the correspondent matched filter. This enables the use of the
radar in critical environments where interception and jamming robustness is a
requisite. Additionally, mutual interference between two radars that occupy
the same transmit spectral band can be made negligible K. Savci et al (2020)
Palo et al. (2020).

It is also possible to implement pulse-to-pulse diversity in ATLAS by loading
a different waveform during the idle time of the previous pulse. Trains of
diverse pulses can lead to reduction of the range sidelobes of the
autocorrelation function. The Golay complementary sequences, for example, can
be used to phase-modulate the carrier wave leading to trains of complementary
pulses. The sum of the autocorrelation of complementary pulses is zero except
for the zero shift, leading to an impulse like response N.Levanon and
E.Mozeson (2004).

Since the system can be used in passive bistatic mode, it will also be used to
test techniques for improved detection by illumination matching on receive
such as the one published in Marques (2018).

##  5 Conclusions and Future Work

In this work we started by doing a review on the current tracking radar
systems deployed around the world as well as the networks responsible for SSA
and SST tasks (Section 1). Next, we explained the necessity of developing a
tracking radar in Portugal and presented the current one in development,
ATLAS. An extensive technical description of the system was provided (Section
2) and simulation studies on its expected performance were developed and
discussed (Section 3). Since waveform design plays a big role in radar
systems, we gave a description of some interesting waveforms that can be used
by the system in order to increase the tracking quality (Section 4).

The system is expected to track targets up to 1000 km in distance with a
minimum RCS of 10 cm2, with an average of 500-700 objects passing hourly in
the visible sky. With proper waveform design and processing, it is possible to
achieve measurements at 1000 km of range with less than 0.05% error.

The ATLAS pilot radar was developed with cutting edge hardware technology and
uses a highly modular architecture with processing fully in the digital
domain. This enables a cost reduction in the development and maintenance of
the system and provides a platform for research and innovation in the radar
field. We also indicated several emerging fields that can be tested and
developed using it, such as matched illumination, noise radar and cognitive
radar (Section 4.2).

Since the system is at an early stage, several tests and upgrades are
envisioned:

  * •

Calibration campaigns for the antenna motorized tracking system;

  * •

Operational testing in real scenarios;

  * •

Incorporation of the system into an operational SST network;

  * •

Addition of a waveform generator with In-Phase and Quadrature Modulation;

conceptualization, João Pandeirada, Miguel Bergano, Domingos Barbosa and Paulo
Marques; investigation and writing, João Pandeirada, Miguel Bergano, João
Neves and Paulo Marques; simulations, results and visualization, João
Pandeirada, Miguel Bergano, Bruno Coelho and Valério Ribeiro; project
administration, Miguel Bergano and Domingos Barbosa; funding acquisition,
Domingos Barbosa and Valério Ribeiro.

The team acknowledges financial support from the Aga Khan Development Network
and the Fundação para a Ciência e a Technologia, Portugal, for the Science and
Technology Cooperation DOPPLER - DevelOpment of PaloP knowLEdge in
Radioastronomy, project number 333197717 . V.A.R.M.R. acknowledges financial
support from the Fundacão para a Ciência e Tecnologia (FCT) in the form of an
exploratory project of reference IF/00498/2015 and PHOBOS, POCI-01-0145-
FEDER-029932, funded by Programa Operacional Competitividade e
Internacionalizacão (COMPETE 2020) and FCT, Portugal. The team acknowledges
financial support from Enabling Green E-science for the Square Kilometre Array
Research Infrastructure (ENGAGE-SKA), POCI-01-0145-FEDER-022217, funded by
COMPETE 2020 and FCT, Portugal and from the European Commission H2020
Programme under the grant agreement 2-3SST2018-20.

###### Acknowledgements.

The team acknowledges the European Commission Horizon 2020 Programme under the
Space Surveillance and Tracking (SST) grant agreement 2-3SST2018-20. The
authors declare no conflict of interest. References

## References

  * (1) Statista.  Number of satellites in orbit by major country as of March 31, 2020.  https://www.statista.com/statistics/264472/number-of-satellites-in-orbit-by-operating-country/#statisticContainer.  [Accessed January, 8, 2021].
  * (2) United Nations Office for Outer Space Affairs.  Treaty on Principles Governing the Activities of States in the Exploration and Use of Outer Space, including the Moon and Other Celestial Bodies.  https://www.unoosa.org/oosa/en/ourwork/spacelaw/treaties/introouterspacetreaty.html.  [Accessed January, 8, 2021].
  * Skolnik (2001) Skolnik, M.L.  Introduction to RADAR Systems; Tata McGraw-Hill, 2001.
  * Chen (2004) Chen, W.K.  The Electrical Engineering Handbook; Elsevier Academic Press, 2004\.
  * Galati (2016) Galati, G.  100 Years of Radar; Switzerland: Springer, 2016.
  * Ingwersen and Lemnios (2000) Ingwersen, P.; Lemnios, W.  Radars for Ballistic Missile Defense Research,.  LINCOLN LABORATORY JOURNAL 2000, 2, 22.
  * Schrogl et al. (2015) Schrogl, K.U.; Hays, P.; Robinson, J.; Moura, D.; Gianopapa, C.  Handbook of Space Security; New York: Springer, 2015.
  * (8) Thales Group.  Stir - Tracking and illumination radar.  https://www.thalesgroup.com/en/stir-tracking-and-illumination-radar.  [Accessed September, 17, 2020].
  * (9) SKA.  SQUARE KILOMETRE ARRAY.  https://www.skatelescope.org/.  [Accessed September, 18, 2020].
  * (10) eoPortal.  Space Fence for SSA (Space Situational Awareness) Services.  https://directory.eoportal.org/web/eoportal/satellite-missions/s/space-fence.  [Accessed September, 18, 2020].
  * Stokely et al. (2006) Stokely, C.L.; Foster, J.L.; Stansbery, E.G.; Benbrook, J.R.; Juarez, Q.  Haystack and HAX Radar Measurements of the Orbital Debris Environment 2006.  p. 242.
  * Murray et al. (2019) Murray, J.; Miller, R.; Matney, M.; Kennedy, T.  Orbital Debris Radar Measurements from the Haystack Ultra-wideband Satellite Imaging Radar (HUSIR): 2014-2017.  First International Orbital Debris Conference 2019.
  * Du (2017) Du, R.  China’s approach to space sustainability: Legal and policy analysis.  Space Policy 2017, 42, 8–16.
  * Prasad (2005) Prasad, M.Y.S.  Technical and legal issues surrounding space debris— India’s position in the UN.  Space Policy 2005, 21, 243–249.
  * Adimurthy and Ganeshan (2006) Adimurthy, V.; Ganeshan, A.  Space debris mitigation measures in India.  Acta Astronautica 2006, 58, 168 – 174.  doi:blackhttps://doi.org/10.1016/j.actaastro.2005.09.002.
  * Froehlich et al. (2020) Froehlich, A.; Alonso, D.; Soria, A.; Marchi, E.D.  Space Supporting Latin America, Latin America’s Emerging Space Middle Power; Springer, 2020; pp. 197–273.  doi:black10.1007/978-3-030-38520-0.
  * Nogueira et al. (2012) Nogueira.; Andrei, E.; Neto, A.; Tang, D.; Li, Z.; Mao, Y.; Penna, Y.; Teixeira, J.; Yu, R.; Neto, Y.; Messias.  The China-Brazil Program of Space Debris Monitoring.  39th COSPAR Scientific Assembly 2012, E1.6-2-12, 1381.
  * DEFENSE INTELLIGENCE AGENCY (2019) (DIA) DEFENSE INTELLIGENCE AGENCY (DIA).  Challenges to Security in Space; 2019; p. 46.
  * Mehrholz (1997) Mehrholz, D.  Radar observations in low earth orbit.  Advances in Space Research 1997, 19, 203–212.  doi:black10.1016/S0273-1177(97)00002-1.
  * Mehrholz et al. (1999) Mehrholz, D.; Leushacke, L.; Jehn, R.  COBEAM-1796 Experiment.  Advances in Space Research 1999, 23, 23–32.
  * Muntoni et al. (2017) Muntoni, G.; Schirru, L.; Pisanu, T.; Montisci, G.; Valente, G.; Gaudiomonte, F.; Serra, G.; Urru, E.; Ortu, P.; Fanti, A.  Space Debris Detection in Low Earth Orbit with the Sardinia Radio Telescope.  Electronics 2017, 6.  doi:black10.3390/electronics6030059.
  * European Union (2014) European Union.  DECISION No 541/2014/EU OF THE EUROPEAN PARLIAMENT AND OF THE COUNCIL of 16 April 2014 establishing a Framework for Space Surveillance and Tracking Support, 2014.
  * (23) EU SST.  EUSST-EUROPEAN SPACE SURVEILLANCE AND TRACKING PROJECTS.  https://www.eusst.eu.
  * Eshbaugh et al. (2014) Eshbaugh, J.V.; Morrison, R.L.; Hoen, E.W.; Hiett, T.C.; Benitz, G.R.  HUSIR Signal Processing.  Lincoln Laboratory Journal 2014, 21, 115–134.
  * Rejto (2000) Rejto, S.  Radar open systems architecture and applications.  Record of the IEEE 2000 International Radar Conference [Cat. No. 00CH37037], 2000, pp. 654–659.  doi:black10.1109/RADAR.2000.851911.
  * Fonseca et al. (2006) Fonseca, R.; Barbosa, D.; Cupido, L.; Mourao, A.; Santos, D.; Smoot, G.; Tello, C.  Site evaluation and RFI spectrum measurements in Portugal at the481frequency range 0.408–10 GHz for a GEM polarized galactic radio emission experiment,.  New Astronomy - NEW ASTRON 2006, 11, 551–556.
  * Cree () (2018) Cree (2018).  CGHV59350 Data Sheet, 350 W, 5200-5900 MHz GaN HEMT for C-Band Radar.  https://www.wolfspeed.com/downloads/dl/file/id/463/product/174/cghv59350.pdf.  [accessed May 29, 2018].
  * Simons et al. (2019) Simons, R.N.; Wintucky, E.G.; Waldstein, S.W.  A Novel Reconfigurable GaN Based Fully Solid-State Microwave Power Module for Communications/Radar Applications.  2019 IEEE Aerospace Conference, 2019, pp. 1–7.  doi:black10.1109/AERO.2019.8741973.
  * (29) General Microwave.  Series D195 Octave-Band PIN Diode Attenuator/Modulators.
  * (30) Wolfspeed Cree.  CGHV59070 70 W, 4.4 - 5.9 GHz, 50 V, RF Power GaN HEMT.  [Rev 2.2 - April 2020].
  * (31) MTI-Milliren Technologies, Inc..  260 Series OCXO.  [Mar 2015].
  * (32) Advantech.  MIO-2360 IntelR©Pentium N4200/Celeron N3350/AtomTME3900 series Pico-ITX SBC.  [Last updated: 14-Sep-2018].
  * Skolnik (2008) Skolnik, M.L.  "An Introduction and Overview of Radar", in "Radar Handbook"; McGraw-Hill, 2008.
  * (34) International Telecommunication Union.  "Propagation data and prediction methods required for the design of Earth-space telecommunication systems", in "Recommendation ITU-R".
  * Riché et al. (2014) Riché, V.; Méric, S.; Baudais, J.; Pottier, î  Investigations on OFDM Signal for Range Ambiguity Suppression in SAR Configuration.  IEEE Transactions on Geoscience and Remote Sensing 2014, 52, 4194–4197.  doi:black10.1109/TGRS.2013.2280190.
  * (36) USSPACECOM.  Space-Track.  https://www.space-track.org/.
  * (37) Satellogic SA.  Orbit Predictor.  https://github.com/satellogic/orbit-predictor.  [Accessed December, 30, 2020].
  * N.Levanon and E.Mozeson (2004a) N.Levanon.; E.Mozeson.  "Ambiguity Function", in "Radar Signals"; John Wiley & Sons, Inc., 2004; pp. 34–52.
  * N.Levanon and E.Mozeson (2004b) N.Levanon.; E.Mozeson.  "Matched Filter", in "Radar Signals"; John Wiley & Sons, Inc., 2004; pp. 20–33.
  * Skolnik (2008) Skolnik, M.L.  "Pulse Compression Radar", in "Radar Handbook"; McGraw-Hill, 2008.
  * N.Levanon and E.Mozeson (2004) N.Levanon.; E.Mozeson.  "Phase-Coded Pulse", in "Radar Signals"; John Wiley & Sons, Inc., 2004; pp. 100–167.
  * Dzvonkovskaya and Rohling (2008) Dzvonkovskaya, A.; Rohling, H.  Long binary phase codes with good autocorrelation properties.  2008 International Radar Symposium, 2008, pp. 1–4.  doi:black10.1109/IRS.2008.4585750.
  * N.Levanon and E.Mozeson (2004) N.Levanon.; E.Mozeson.  "Basic Radar Signals", in "Radar Signals"; John Wiley & Sons, Inc., 2004; pp. 53–73.
  * Marques (2017) Marques, P.A.  Noise radar detection optimized for selected targets.  2017 18th International Radar Symposium (IRS), 2017, pp. 1–9.  doi:black10.23919/IRS.2017.8008226.
  * K. Savci et al (2020) K. Savci et al.  Noise Radar—Overview and Recent Development.  IEEE Aerospace and Electronic498Systems Magazine 2020, 35, 8–20.  doi:black10.1109/MAES.2020.2990591.
  * Palo et al. (2020) Palo, F.D.; Galati, G.; Pavan, G.; Wasserzier, C.; Savci, K.  Introduction to Noise Radar and Its Waveforms.  Sensors 2020, 20.  doi:black10.3390/s20185187.
  * N.Levanon and E.Mozeson (2004) N.Levanon.; E.Mozeson.  "Coherent train of diverse pulses", in "Radar Signals"; John Wiley & Sons, Inc., 2004; pp. 262–270.
  * Marques (2018) Marques, P.A.  More Illumination for Passive Radar.  European Conf. on Synthetic Aperture Radar, 2018, Vol. 1.