#  Skip-sampling: subsampling in the frequency domain

Tucker S. McElroy
U.S. Census Bureau
4600 Silver Hill Road, Washington, DC 20233
<EMAIL_ADDRESS>&Dimitris N. Politis
Department of Mathematics and Halicioglu Data Science Institute
University of California at San Diego
La Jolla, CA 92093-0112, USA<EMAIL_ADDRESS>
###### Abstract

Over the last 35 years, several bootstrap methods for time series have been
proposed. Popular ‘time-domain’ methods include the block-bootstrap, the
stationary bootstrap, the linear process bootstrap, etc.; subsampling for time
series is also available, and is closely related to the block-bootstrap.
‘Frequency-domain’ bootstrap has been performed either by resampling the
periodogram ordinates or by resampling the ordinates of the Discrete Fourier
Transform (DFT). The paper at hand proposes a novel construction of
subsampling the DFT ordinates, and investigates its theoretical properties and
realm of applicability.

_Keywords_ Discrete Fourier Transform, Spectral Density, Time Series

##  1 Introduction

Efron (1979) developed the bootstrap for independent and identically
distributed (i.i.d.) data, and paved the way for practical nonparametric
statistics in the modern era. Soon after, practitioners were able to apply
resampling ideas in a variety of non-i.i.d. situations including the
interesting case of dependent data.

Time series analysis has a ‘time-domain’ and a ‘frequency-domain’ aspect to
it. Consequently, proposals for resampling time series can be in either of
these two flavors. There have been several proposals with regards to ‘time-
domain’ resampling plans; these include the block-bootstrap in its many
variations, the stationary bootstrap, the linear process bootstrap, etc.; see
Chapter 12 of McElroy and Politis (2020) for a description. Subsampling for
time series is closely related to the block-bootstrap; see Politis et al.
(1999). For book-length treatments of the state-of-the-art in resampling time
series, see Lahiri (2003a) or Kreiss and Paparoditis (2023).

One of the first papers on a ‘frequency-domain’ bootstrap was by Franke and
Härdle (1992), who proposed resampling the periodogram ordinates. The
motivation behind that approach is that periodogram ordinates at different
Fourier frequencies are approximately independent. Different avenues based on
this approach were pursued by several researchers including Janas and Dahlhaus
(1994), Dahlhaus and Janas (1996), Kreiss and Paparoditis (2003), and Meyer et
al. (2020).

Interestingly, there was an earlier—albeit unpublished—report by Hurvich and
Zeger (1987), who proposed resampling the ordinates of the Discrete Fourier
Transform (DFT), as they are also approximately independent. Actually, the
aforementioned (approximate) independence of periodogram ordinates is a
consequence of the (approximate) independence of DFT ordinates, since the
periodogram is a function of the DFT. Hence, resampling the DFT can be thought
of as a more fundamental construction; a rigorous development can be found in
Kirch and Politis (2011).

Since resampling the DFT is a fundamental construction, the question of
possibly subsampling the DFT presents itself; this is the subject of the paper
at hand. The basic idea is to divide the DFT (based on a sample of size $T$)
into $q$ vectors of length $b$, each consisting of the DFT ordinates at
frequencies separated by $q/T$. Each such vector is asymptotically independent
of one another, and distributed as a DFT vector based on a sample of size $b$.
If the statistic at hand is computable based on the DFT alone, it could be re-
computed on the smaller DFTs, and an empirical distribution of such subsample
statistics, appropriately centered and normalized, would estimate the original
statistic’s sampling distribution.

The above construction will be termed skip-sampling of the DFT because of the
process of skipping over some frequencies in putting together the subsample
DFT vectors. The following section gives the precise construction as well as
some background on the properties of the DFT. Some theoretical results on
skip-sampling of the DFT are given in Section 3, while applications to
spectral means and ratio statistics are given in Section 4.

##  2 Problem setup

Let $X_{1},X_{2},\ldots,X_{T}$ be an observed sample from a strictly
stationary time series $\\{X_{t}\\}$, and denote ${\bf
X}={[X_{1},X_{2},\ldots,X_{T}]}^{\prime}$. Time series analysis in the
‘frequency domain’ hinges on the Discrete Fourier Transform (DFT), which maps
the data vector ${\bf X}$ to a vector with (approximately) independent
entries. To define the DFT, consider the Fourier frequencies
$\lambda_{\ell}=2\pi\ell/T$, where $\ell$ is an integer satisfying
$[T/2]-T+1\leq\ell\leq[T/2]$ and $[\cdot]$ denotes integer part; this index
range corresponds to $-[T/2]\leq\ell\leq[T/2]$ when $T$ is odd, or
$-[T/2]+1\leq\ell\leq[T/2]$ when $T$ is even.111Some authors define the DFT
using the Fourier frequencies $\lambda_{\ell}=2\pi\ell/T$ for
$\ell=0,1,\ldots,T-1$; due to the periodicity present in the DFT, there is no
discrepancy here, other than a re-ordering of the entries of the DFT vector
which will help us exploit some symmetry properties. Define a $T\times T$
matrix $Q$ with complex-valued entries
$Q_{jk}=T^{-1/2}\,e^{ij\lambda_{[T/2]-T+k}}$ for $1\leq j,k\leq T$. Note that
$Q$ is unitary, i.e., $Q^{-1}=Q^{*}$, the conjugate transpose. The DFT vector
(see Proposition of 7.2.7 of McElroy and Politis (2020) for more details) is
defined as $\widetilde{{\bf X}}=Q^{*}\,{\bf X}$, which means that the $k$th
component of the DFT vector is

$\widetilde{X}_{j}=T^{-1/2}\sum_{k=1}^{T}e^{-ik\lambda_{[T/2]-T+j}}X_{k}.$ (1)

The DFT map is invertible, because clearly ${\bf X}=Q\,\widetilde{{\bf X}}$.

###  2.1 DFT Symmetries

The DFT vector has certain symmetries; in order to describe these symmetries,
we define a transposition matrix $P$ (or $P_{T}$, when we need to annotate its
dimension), which is a $T\times T$ matrix with ones on the trans-diagonal (and
zeros elsewhere); its action on a vector is to reverse the order of its
components. Let $1_{T}$ denote the identity matrix of dimension $T$, and let
$\Pi$ be the permutation matrix that when applied to a column vector shifts
all the components upwards one position, and sends the first component to the
last (bottom) position. The complex conjugate of $z$ is denoted
$\mathcal{C}z$, and the real and imaginary parts are $\Re z$ and $\Im z$
respectively.

###### Proposition 1

If $T$ is odd, the DFT vector $\widetilde{{\bf X}}$ satisfies

$P\,\widetilde{\bf X}=\mathcal{C}{\widetilde{\bf X}}.$ (2)

If $T$ is even, the DFT vector satisfies

$\Pi\,P\,\widetilde{\bf X}=\mathcal{C}{\widetilde{\bf X}}.$ (3)

#### Proof of Proposition 1.

First consider the case that $T$ is odd, so $T=2m+1$ for some integer $m$.
Then $P\ \widetilde{\bf X}$ has entries in reverse order, so that the middle
component is unchanged but all others are flipped. Using (1),
$\widetilde{{X}}_{j}=T^{-1/2}\sum_{k=1}^{T}e^{-ik\lambda_{-m-1+j}}X_{k}$. On
the other hand, the $j$th component of $P\ \widetilde{\bf X}$ is
$\widetilde{{X}}_{T+1-j}=T^{-1/2}\sum_{k=1}^{T}e^{-ik\lambda_{m+1-j}}X_{k}$,
which is the conjugate of $\widetilde{X}_{j}$. This proves (2). Next, suppose
that $T$ is even, so $T=2m$ for some integer $m$. For $1\leq j\leq T-1$, the
$j$th component of $\Pi\,P\,\widetilde{\bf X}$ is the $j+1$th component of
$P\,\widetilde{\bf X}$, which is

$\widetilde{X}_{T-j}=T^{-1/2}\sum_{k=1}^{T}e^{-ik\lambda_{m-j}}X_{k}=T^{-1/2}\sum_{k=1}^{T}e^{ik\lambda_{-m+j}}X_{k}=\mathcal{C}\widetilde{X}_{j}.$

Moreover, the $T$th component of $\Pi\,P\,\widetilde{\bf X}$ is the first
component of $P\,\widetilde{\bf X}$, i.e., $\widetilde{X}_{T}$. Because
$\lambda_{m}=2\pi m/T=\pi$, this number is real, and hence
$\widetilde{X}_{T}=\mathcal{C}\widetilde{X}_{T}$. This proves (3). $\quad\Box$

In view of Proposition 1, we say that the DFT vector satisfies a “Symmetry
Property,” defined as follows.

###### Definition 1

A length $T$ complex vector $\widetilde{\bf X}$ satisfies the Symmetry
Property if and only if (2) holds when $T$ is odd and (3) holds when $T$ is
even.

A general length $T$ complex vector ${\bf Z}$ may satisfy the Symmetry
Property, in which case necessarily the vector has a particular structure. If
${\bf Z}$ satisfies (2), it must be the case that

$\Re[P{\bf Z}]=\Re[{\bf Z}],\qquad\Im[P{\bf Z}]=-\Im[{\bf Z}].$

Because the middle entry of an odd-length ${\bf Z}$ has its value unchanged
after application of $P$, its value must be real. Similarly, a vector ${\bf
Z}$ satisfying (3) has real entries for components $T/2$ and $T$, whereas the
subvector of components $1$ through $T/2-1$ is the complex conjugate of the
transposition of the subvector for components $T/2+1$ through $T-1$. In this
paper we will be constructing DFT vectors in various ways, but we need to
ensure that these constructions have the correct properties. In particular,
given a complex vector ${\bf Z}$ it behooves us to know how it can be modified
so as to have the Symmetry Property. The following proposition justifies this
motivation.

###### Proposition 2

If ${\bf Z}$ has the Symmetry Property, then $Q{\bf Z}$ has real-valued
entries.

#### Proof of Proposition 2.

First, noting that $P_{jk}=1_{\\{j+k=T+1\\}}$, we obtain

$\displaystyle{\\{\mathcal{C}QP\\}}_{jk}$
$\displaystyle=\sum_{\ell=1}^{T}\mathcal{C}{Q}_{j\ell}\,P_{\ell
k}=\overline{Q}_{j,T+1-k}$ $\displaystyle=T^{-1/2}\,\exp\\{-2\pi
ij([T/2]-T+T+1-k)/T\\}$ $\displaystyle=T^{-1/2}\,\exp\\{2\pi
ij(k-1-[T/2])/T\\}.$

If $T$ is odd, then $-1-[T/2]=[T/2]-T$, and
${\\{\mathcal{C}{Q}P\\}}_{jk}=Q_{jk}$. But if $T$ is even, then
$-[T/2]=[T/2]-T$ and ${\\{\mathcal{C}{Q}P\\}}_{jk}=Q_{j,k-1}$ for $2\leq k\leq
T$; also ${\\{\mathcal{C}{Q}P\\}}_{j1}=Q_{j,T}$ because $\exp\\{2\pi
ij([T/2]-T)/T\\}=\exp\\{2\pi ij([T/2])/T\\}$. Therefore, when $T$ is odd
$\mathcal{C}{Q}P=Q$, but when $T$ is even $\mathcal{C}{Q}P=Q\Pi$. Next,
because $P$ is idempotent

$\mathcal{C}[Q\,{\bf Z}]=\mathcal{C}[Q]P\,P\mathcal{C}[{\bf
Z}]=\begin{cases}Q\,P\mathcal{C}[{\bf Z}]\quad\mbox{if}\;T\;\mbox{is odd}\\\
Q\,\Pi P\mathcal{C}[{\bf Z}]\quad\mbox{if}\;T\;\mbox{is even}\end{cases}=Q{\bf
Z},$

using (2) and (3). Hence $Q{\bf Z}$ equals its own conjugate, and therefore
must be real. $\quad\Box$

###### Remark 2.1

As an application, we can alter a given complex vector ${\bf Z}$ to have the
Symmetry Property as follows. If the length $T$ is odd, replace the first
$[T/2]$ entries with the conjugate of the first $[T/2]$ entries of $P{\bf Z}$,
and discard the imaginary part of the middle entry in position $[T/2]+1$. If
the length of $T$ is even, we replace components $1$ through $T/2-1$ with the
conjugate of components $T-1$ through $T/2+1$ (so their order is flipped);
also, the imaginary portions of components $T/2$ and $T$ are discarded. These
operations ensure that the modified ${\bf Z}$ has the Symmetry Property.

###  2.2 Subsampling the DFT

We now provide details on the novel construction that is at the heart of this
paper’s methodology. For simplicity, consider positive integers $q$ and $b$
such that222In practical applications, if such an exact choice is not
feasible, then one could let $q=[T/b]$, and work as if the data were just
$X_{1},\ldots,X_{bq}$.  $b\,q=T$, and define sub-components of the DFT vector
by

$\widetilde{{\bf
X}}^{(j)}={[\widetilde{{X}}_{j},\widetilde{{X}}_{q+j},\ldots,\widetilde{{X}}_{(b-1)q+j}]}^{\prime}$

for $j=1,\ldots,q$. In terms of the entire DFT, this operation can be
expressed as $\widetilde{{\bf X}}^{(j)}=[1_{b}\otimes
e_{j}^{\prime}]\,\widetilde{{\bf X}}$, with $e_{j}$ the $j$th unit vector in
$\mathbb{R}^{q}$ and $\otimes$ is the Kronecker product. Because of the
construction of keeping every $q$th Fourier frequency and skipping over the
intervening ones, this operation can be called skip-sampling on the DFT, and
$\widetilde{{\bf X}}^{(j)}$ is called the $j$th skip-sample DFT; it is a
complex vector of length $b$, obtained by evaluating the DFT only at
frequencies of the form $2\pi(\ell q+j)/T$, where
$\ell=[b/2]-b+1,\ldots,[b/2]$.

Recall that the DFT $\widetilde{{\bf X}}$ contains all the information carried
in the sample ${\bf X},$ since we can re-create ${\bf X}$ as
$Q\,\widetilde{{\bf X}}$. However, the $j$th skip-sample DFT $\widetilde{{\bf
X}}^{(j)}$ contains only a part of the information carried by the sample ${\bf
X}$; putting all the skip-sample DFTs $\widetilde{{\bf X}}^{(j)}$ together for
$j=1,\ldots,q$, we can capture the whole information again. In this sense,
working with the skip-sample DFTs $\widetilde{{\bf X}}^{(j)}$ for
$j=1,\ldots,q$ can be considered a form of subsampling in the frequency
domain; this should be contrasted to the usual subsampling of a time series in
the time domain which is done by carving the sample $X_{1},\ldots,X_{T}$ into
smaller blocks, each consisting of $b$ consecutive data points—see Politis et
al. (1999).

Note that $\widetilde{{\bf X}}^{(j)}$ will not necessarily have the Symmetry
Property; employing the techniques of Remark 2.1, we ensure that applying the
$b$-dimensional version of matrix $Q$ to the symmetrized skip-sample DFT (so
as to invert the DFT and bring us back to the time domain) will yield a real-
valued vector of length $b$; this can be useful for statistics that are
formulated in the time domain. However, there is an interesting class of
statistics that are defined in the frequency domain; three prime examples are
discussed in Section 4. The next section defines the new skip-sampling
methodology for such statistics defined in the frequency domain.

##  3 Skip-sampling: the basic methodology and some key results

As before, let $X_{1},X_{2},\ldots,X_{T}$ be an observed sample from a
strictly stationary time series $\\{X_{t}\\}$ with mean $\mu$ and absolutely
summable autocovariance $\gamma_{k}=\mbox{Cov}(X_{0},X_{k})$ so that the
spectral density $f(\lambda)=\sum_{-\infty}^{\infty}\gamma_{k}e^{-ik\lambda}$
is well-defined and continuous on $[-\pi,\pi]$, i.e., belongs to
$C[-\pi,\pi]$. A crude estimate of $f(\lambda)$ is given by the periodogram
$I_{T}(\lambda)=\sum_{-T+1}^{T-1}\hat{\gamma}_{k}e^{-ik\lambda}$, where the
sample autocovariance is defined as
$\hat{\gamma}_{k}=T^{-1}\sum_{t=1}^{T-|k|}(X_{t}-\bar{X})(X_{t+|k|}-\bar{X})$,
and $\bar{X}=T^{-1}\sum_{t=1}^{T}X_{t}$ is the sample mean.

Interestingly, when evaluated at a (nonzero) Fourier frequency, the
periodogram equals the squared magnitude of the DFT. To see that, note the
identity
$I_{T}(\lambda)=T^{-1}|\sum_{t=1}^{T}(X_{t}-\bar{X})e^{-it\lambda}|^{2}$. One
of the columns of the matrix $Q$ consists of constant elements, and the other
columns are orthogonal to it. Hence,

$I_{T}(0)=0,\ \ \mbox{and when}\ \ \ell\neq 0\ \ \mbox{we have}\ \
I_{T}(\lambda_{\ell})=T^{-1}|\sum_{t=1}^{T}X_{t}e^{-it\lambda_{\ell}}|^{2},\ \
\mbox{where}\ \ \lambda_{\ell}=2\pi\ell/T.$

Let $\theta$ be a parameter of interest; we will assume that $\theta$ is some
real-valued functional of $f$, i.e., that $\theta={\cal G}(f)$ for some ${\cal
G}:C[-\pi,\pi]\to{\bf R}$. The periodogram is asymptotically unbiased but
inconsistent for $f(\lambda)$, as its variance does not tend to zero; see
Chapter 9 of McElroy and Politis (2020). However, there are several situations
where $\theta$ can be consistently estimated using the periodogram as a basis.
So consider a statistic $\widehat{\theta}_{T}$ that is a functional of
$I_{T}$, i.e., that $\widehat{\theta}_{T}={\cal G}_{T}(I_{T})$ where, for each
$T$, we have ${\cal G}_{T}:C[-\pi,\pi]\to{\bf R}$. In simpler cases, the
functional ${\cal G}_{T}$ might not depend on $T$, as in the case of spectral
means and ratio statistics discussed in Sections 4.1 and 4.2.

In terms of feasible statistical computing, we will further assume —as it is
invariably the case— that the statistic $\widehat{\theta}_{T}$ is computable
based on the periodogram evaluated just on the Fourier frequencies. Since the
periodogram evaluated at (nonzero) Fourier frequencies equals the squared
magnitude of the DFT, we will assume that

$\widehat{\theta}_{T}={{\cal H}}_{T}\left(\widetilde{{\bf X}}\right),$ (4)

where for each $T$, the function ${{\cal H}}_{T}$ maps ${\bf C}^{T}$ to ${\bf
R}$. We will further assume:

Assumption (A): For some nondegenerate limit distribution $J$, we have
$a_{T}(\widehat{\theta}_{T}-\theta)\stackrel{{\scriptstyle{\cal
L}}}{{\Longrightarrow}}J$ as $T\to\infty$, where $a_{T}=T^{\delta}L(T)$ for
some $\delta>0$ and some slowly varying function $L$.

Letting $J_{T}(x)=Prob\\{a_{T}(\widehat{\theta}_{T}-\theta)\leq x\\}$,
Assumption (A) implies that $J_{T}(x)\to J(x)$ for all points $x$ at which $J$
is continuous.

We can now define the $j$th skip-sample statistic

$\widehat{\theta}_{b}^{(j)}={{\cal H}}_{b}\left(\widetilde{{\bf
X}}^{(j)}\right)$ (5)

for $j=1,\ldots,q$. The idea is that $\widehat{\theta}_{b}^{(j)}$ will have
the same asymptotic distribution as $\widehat{\theta}_{b}^{(1)}$ when
$b\to\infty$. Futhermore, under standard conditions —see e.g. Lahiri (2003b)—,
the DFT ordinates evaluated at different Fourier frequencies will be
asymptotically independent; this would render the skip-sample statistics
$\widehat{\theta}_{b}^{(1)},\ldots,\widehat{\theta}_{b}^{(q)}$ (for fixed $q$)
approximately independent as well.

We formulate these stylized facts in the following Assumption, which operates
under the condition

$\frac{b}{T}+\frac{1}{b}\to 0\ \ \mbox{as}\ \ T\to\infty.$ (6)

Assumption (A∗): Under condition (6) the following are true: (a) For any $j$,
$Prob\\{a_{b}(\widehat{\theta}_{b}^{(j)}-\theta)\leq x\\}-J_{b}(x)=o(1)$ for
all points $x$ at which $J$ is continuous; and (b) for any $j\neq k$, and any
bounded functions $g_{1},g_{2}$, we have
$\mbox{Cov}(g_{1}(\widehat{\theta}_{b}^{(j)}),g_{2}(\widehat{\theta}_{b}^{(k}))\to
0$.

The quantity $a_{T}(\widehat{\theta}_{T}-\theta)$ is sometimes called a
‘root’. Our core result is a “subsampling in the frequency domain” consistency
theorem that gives conditions under which the empirical distribution of the
skip-sample roots $a_{b}(\widehat{\theta}_{b}^{(j)}-\widehat{\theta}_{T})$ for
$j=1,\ldots,q$ can be used to approximate the distribution of the original
root. To develop it, define the two skip-sampling distributions:

$U_{b,T}(x)=q^{-1}\sum_{j=1}^{q}{\bf
1}\\{a_{b}(\widehat{\theta}_{b}^{(j)}-\theta)\leq x\\}\ \ \mbox{and}\ \
L_{b,T}(x)=q^{-1}\sum_{j=1}^{q}{\bf
1}\\{a_{b}(\widehat{\theta}_{b}^{(j)}-\widehat{\theta}_{T})\leq x\\},$

where ${\bf 1}$ denotes the indicator function. Of the two skip-sampling
distributions, $U_{b,T}$ is termed an ‘oracle’ as it requires knowledge of
$\theta$ for its construction. By contrast, $L_{b,T}$ is a bona fide statistic
that can be used for estimation purposes.

###### Theorem 1

Assume condition (6) and Assumptions (A) and (A∗). Then,

$L_{b,T}(x)\stackrel{{\scriptstyle P}}{{\longrightarrow}}J(x)$

for all points $x$ at which $J$ is continuous; here $\stackrel{{\scriptstyle
P}}{{\longrightarrow}}$ denotes convergence in probability.

Proof. Let $x$ be a point of continuity of $J$. By the first argument given in
the proof of Theorem 2.2.1 of Politis et al. (1999), $U_{b,T}(x)$ and
$L_{b,T}(x)$ are asymptotically close. So to prove Theorem 1, it suffices to
show that $U_{b,T}(x)\stackrel{{\scriptstyle P}}{{\longrightarrow}}J(x)$.

Note that
$EU_{b,T}(x)=q^{-1}\sum_{j=1}^{q}Prob\\{a_{b}(\widehat{\theta}_{b}^{(j)}-\theta)\leq
x\\}=J_{b}(x)+o(1)$ by Assumption (A∗) and the Cesaro sums lemma. Furthermore,
$J_{b}(x)\to J(x)$ by Assumption (A) and (6). Now note that

$\mbox{Var}(U_{b,T}(x))=q^{-2}\sum_{j=1}^{q}\sum_{k=1}^{q}\mbox{Cov}(g_{1}(\widehat{\theta}_{b}^{(j)}),g_{1}(\widehat{\theta}_{b}^{(k})),$

where $g_{1}(\widehat{\theta}_{b}^{(j)})={\bf
1}\\{a_{b}(\widehat{\theta}_{b}^{(j)}-\theta)\leq x\\}$. By Assumption (A∗)
and Cesaro sums, it follows that

$\mbox{Var}(U_{b,T}(x))=q^{-2}\sum_{j=1}^{q}\mbox{Cov}(g_{1}(\widehat{\theta}_{b}^{(j)}),g_{1}(\widehat{\theta}_{b}^{(j)}))+o(1).$

Since $g_{1}$ is an indicator, it follows that
$|\mbox{Cov}(g_{1}(\widehat{\theta}_{b}^{(j)}),g_{1}(\widehat{\theta}_{b}^{(j)}))|\leq
1$. Hence, $\mbox{Var}(U_{b,T}(x))=o(1)$, and the desired result follows by
Chebyshev’s inequality. $\Box$

###### Remark 3.1

In many situations, the limit law $J$ will be $N(0,v)$. In this case, it may
be of interest to use skip-sampling to estimate the asymptotic variance $v$,
and use the normal tables (instead of the quantiles of the skip-sampling
distribution $L_{b,T}$) in order to construct confidence intervals and tests.
The skip-sampling estimator of $v$ is given by

$\widehat{v}_{b}=\frac{a_{b}^{2}}{q}\sum_{j=1}^{q}(\widehat{\theta}_{b}^{(j)}-\bar{\widehat{\theta}}_{b})^{2},$
(7)

where
$\bar{\widehat{\theta}}_{b}=q^{-1}\sum_{j=1}^{q}\widehat{\theta}_{b}^{(j)}$.
The consistency of $\widehat{v}_{b}$ requires some different conditions that
are outlined in Corollary 1 below. Such conditions can be verified for the two
prominent types of periodogram-based statistics, namely spectral means and
ratio statistics; see Sections 4.1 and 4.2. Additional examples of potential
applicabilty of frequency domain resampling (including skip-sampling) are
given in Corollary 3.1 of Bertail and Dudek (2021) and its related discussion.

###### Corollary 1

Assume Assumption (A) with $\sup_{T}E\widehat{\theta}_{T}^{4}<\infty$, and
$a_{T}^{2}\mbox{Var}[\widehat{\theta}_{T}]\to v>0$ as $T\to\infty$. Let $b$ be
a sequence satisfying

$\frac{b}{T}+\frac{1}{b}+\frac{a_{b}^{2}}{T}\to 0\ \ \mbox{as}\ \ T\to\infty.$
(8)

Also assume that, for any $1\leq i,j\leq q$ and $i\neq j$, the following set
of assumptions holds:

$\begin{cases}&\mathbb{E}[\widehat{\theta}_{b}^{(j)}]=\theta+o(a_{b}^{-1})\\\
&a_{b}^{2}\mbox{Var}[\widehat{\theta}_{b}^{(j)}]=v+o(1)\\\
&\mbox{Cov}[\widehat{\theta}_{b}^{(i)},\widehat{\theta}_{b}^{(j)}]=O(T^{-1}).\end{cases}$
(9)

Further assume that when $i\neq j$,

$\mbox{Cov}\left[a_{b}^{2}(\widehat{\theta}_{b}^{(i)}-{\theta})^{2},\
a_{b}^{2}(\widehat{\theta}_{b}^{(j)}-{\theta})^{2}\right]=o(1).$ (10)

Then, $\widehat{v}_{b}\stackrel{{\scriptstyle P}}{{\longrightarrow}}v$ as
$T\to\infty$.

###### Remark 3.2

The validity of assumption (10) can be motivated by the asymptotic
independence of $\widehat{\theta}_{b}^{(i)}$ and $\widehat{\theta}_{b}^{(j)}$.
Moreover, note that the set of assumptions (9) implies
$a_{b}^{2}E(\widehat{\theta}_{b}^{(i)}-{\theta})^{2}=v+o(1)$. Hence, by
Markov’s inequality,
$a_{b}^{2}(\widehat{\theta}_{b}^{(i)}-{\theta})^{2}=v+o_{P}(1)$, implying that
$\left[a_{b}^{2}(\widehat{\theta}_{b}^{(i)}-{\theta})^{2}-v\right]\left[a_{b}^{2}(\widehat{\theta}_{b}^{(j)}-{\theta})^{2}-v\right]$
$=o_{P}(1)$. Eq. (10) can be be viewed as a stronger version of this result.

#### Proof of Corollary 1.

Consider the ‘oracle’ quantity
$\widetilde{v}_{b}=\frac{a_{b}^{2}}{q}\sum_{j=1}^{q}(\widehat{\theta}_{b}^{(j)}-{\theta})^{2},$
and note that

$\widehat{v}_{b}=\frac{a_{b}^{2}}{q}\sum_{j=1}^{q}(\widehat{\theta}_{b}^{(j)}-\theta+\theta-\bar{\widehat{\theta}}_{b})^{2}=\widetilde{v}_{b}-a_{b}^{2}(\bar{\widehat{\theta}}_{b}-\theta)^{2}.$
(11)

By assumption, $E\bar{\widehat{\theta}}_{b}=\theta+o(a_{b}^{-1})$, and

$q^{2}\mbox{Var}[\bar{\widehat{\theta}}_{b}]=\sum_{i=1}^{q}\sum_{j=1}^{q}\mbox{Cov}[\widehat{\theta}_{b}^{(i)},\widehat{\theta}_{b}^{(j)}]=\sum_{i=1}^{q}\mbox{Var}[\widehat{\theta}_{b}^{(i)}]+\sum_{i=1}^{q}\sum_{j\neq
i}\mbox{Cov}[\widehat{\theta}_{b}^{(i)},\widehat{\theta}_{b}^{(j)}]=qa_{b}^{-2}(v+o(1))+O(\frac{q^{2}}{T})$

so that

$a_{b}^{2}E(\bar{\widehat{\theta}}_{b}-\theta)^{2}=\frac{v+o(1)}{q}+O(\frac{a_{b}^{2}}{T}),$

which tends to zero by (8). Hence, eq. (11) implies that
$\widehat{v}_{b}-\widetilde{v}_{b}\stackrel{{\scriptstyle
P}}{{\longrightarrow}}0$.

To show $\widetilde{v}_{b}\stackrel{{\scriptstyle P}}{{\longrightarrow}}v$,
first note that assumption (9) implies

$E\widetilde{v}_{b}=\frac{a_{b}^{2}}{q}\sum_{j=1}^{q}\left(\mbox{Var}[\widehat{\theta}_{b}^{(j)}]+\mbox{Bias}^{2}[\widehat{\theta}_{b}^{(j)}]\right)=v+o(1),$

i.e., $\mbox{Bias}[\widetilde{v}_{b}]\to 0$. To show
$\mbox{Var}[\widetilde{v}_{b}]\to 0$, note that

$\widetilde{v}_{b}-v=\frac{a_{b}^{2}}{q}\sum_{j=1}^{q}(\widehat{\theta}_{b}^{(j)}-{\theta})^{2}-v=\frac{1}{q}\sum_{j=1}^{q}s_{j},$

letting $s_{j}=a_{b}^{2}(\widehat{\theta}_{b}^{(j)}-{\theta})^{2}-v$. Finally,

$E(\widetilde{v}_{b}-v)^{2}=\frac{1}{q^{2}}\sum_{i=1}^{q}\sum_{j=1}^{q}E(s_{i}s_{j})=\frac{1}{q^{2}}\sum_{i=1}^{q}Es_{i}^{2}+\frac{1}{q^{2}}\sum_{i=1}^{q}\sum_{j\neq
i}E(s_{i}s_{j})=O(1/q)+o(1)$

by (10). By Chebyshev’s inequality, it follows that
$\widetilde{v}_{b}\stackrel{{\scriptstyle P}}{{\longrightarrow}}v$, and
therefore $\widehat{v}_{b}\stackrel{{\scriptstyle P}}{{\longrightarrow}}v$ as
well. $\Box$

##  4 Examples of skip-sampling applicability

###  4.1 Spectral means

Consider a bounded function $g(\lambda)$ of domain $[-\pi,\pi]$ that has
bounded variation, and denote $\langle
g\rangle={(2\pi)}^{-1}\,\int_{-\pi}^{\pi}g(\lambda)\,d\lambda.$ A linear
spectral mean (Dahlhaus, 1985) is a parameter of the form

$\theta=\langle
g\,f\rangle=\frac{1}{2\pi}\,\int_{-\pi}^{\pi}g(\lambda)\,f(\lambda)\,d\lambda.$
(12)

The prime example of a linear spectral mean is the autocovariance at lag $k$,
where $g(\lambda)=e^{ik\lambda}$.

As already mentioned, the periodogram is asymptotically unbiased but
inconsistent for $f(\lambda)$, as its variance does not tend to zero. However,
plugging in $I_{T}$ instead of $f$ in eq. (12) yields a consistent estimator
of the spectral mean, since integration works like summation in terms of
reducing the variance. As a matter of fact, the integral in eq. (12) is
typically approximated by a Riemann sum over the Fourier frequencies.
Consequently, a linear spectral mean $\theta$ satisfying eq. (12) is
practically estimated by

$\widehat{\theta}_{T}=T^{-1}\,\sum_{\ell\in
R_{T}}g(2\pi\ell/T)\,I_{T}(2\pi\ell/T),$

where we define the index range $R_{T}=\\{\ell:[T/2]-T+1\leq\ell\leq[T/2]\\}$.

To prove our next result, we consider stationary non-Gaussian processes that
satisfy autocumulant conditions described in Taniguchi and Kakizawa (2000).
Supposing that all moments exist and the autocumulant functions are defined
via

$\gamma_{h_{1},\ldots,h_{k-1}}=\mbox{cum}\\{X_{t+h_{1}},X_{t+h_{2}},\ldots,X_{t+h_{k-1}},X_{t}\\},$

we will entertain:

Assumption (B): for all $k\geq 2$ and each $j=1,\ldots,k-1$ we have

$\sum_{h_{1}\in\mathbb{Z}}\cdots\sum_{h_{k-1}\in\mathbb{Z}}(1+|h_{j}|)\,|\gamma_{h_{1},\ldots,h_{k-1}}|<\infty.$

For a process satisfying Assumption (B) the tri-spectral density $F$ is well-
defined, and is given by

$F(\omega_{1},\omega_{2},\omega_{3})=\sum_{h_{1}\in\mathbb{Z}}\sum_{h_{2}\in\mathbb{Z}}\sum_{h_{3}\in\mathbb{Z}}\gamma_{h_{1},h_{2},h_{3}}\exp\\{-i(h_{1}\omega_{1}+h_{2}\omega_{2}+h_{3}\omega_{3})\\},$

in which case it can be shown (McElroy and Roy, 2022) that

$\sqrt{T}\,\left(\widehat{\theta}_{T}-\theta\right)\stackrel{{\scriptstyle{\cal
L}}}{{\Longrightarrow}}\mathcal{N}\left(0,\langle
g\,g^{\star}\,f^{2}\rangle+\langle\langle g\,g\,F\rangle\rangle\right)\ \
\mbox{as}\ \ T\to\infty;$ (13)

here, we have used the short-hand
$g^{\star}(\lambda)=g(\lambda)+g^{\sharp}(\lambda),$ where $g^{\sharp}$ is the
reflection of $g$ about the y-axis, and we have denoted

$\langle\langle
g\,g\,F\rangle\rangle=\frac{1}{{(2\pi)}^{2}}\,\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}g(\lambda)\,g(\omega)\,F(\lambda,-\lambda,\omega)\,d\lambda
d\omega.$

Dahlhaus (1985) proved eq. (13) for linear processes. Working under Assumption
(B) allows us to go beyond the setting of linearity; however, there could be
different sufficient conditions for (13).

Evaluating the linear spectral mean on the $j$th skip-sample DFT yields

$b^{-1}\,\sum_{\ell\in R_{b}}g(2\pi\ell/b)\,I_{T}(2\pi\ell/b+2\pi j/T).$ (14)

The displacement of the periodogram by $2\pi j/T$ means that this is no longer
an even function of $2\pi\ell/b$ (unless $j=q$), so the asymptotic variance
has no contribution from $g^{\sharp}$. In order to correct this, we need to
define skip-samples over $[0,\pi]$ and then reflect onto $[-\pi,0]$, which
suggests the definition of the $j$th skip-sample statistic

$\widehat{\theta}_{b}^{(j)}=b^{-1}\,\sum_{\ell=1}^{[b/2]}g^{\star}(2\pi\ell/b)\,I_{T}(2\pi\ell/b+2\pi
j/T).$

The above is obtained from (14) by imposing that at Fourier frequencies
$2\pi\ell/b$ with negative $\ell$, we evaluate the periodogram at
$2\pi\ell/b-2\pi j/T$ instead of at $2\pi\ell/b+2\pi j/T$. Then, ignoring the
contribution from $\ell=0$, we obtain the above expression for
$\widehat{\theta}_{b}^{(j)}$ in terms of $g^{\star}(\lambda)$. Next, we
summarize some of the moment properties of these skip-sample statistics.

###### Theorem 2

Assume that $\\{X_{t}\\}$ is strictly stationary and satisfies Assumption B.
Consider a linear spectral mean $\theta$ satisfying eq. (12), and some fixed
bounded function $g(\lambda)$ having bounded variation. Let $b$ be a sequence
satisfying (6). Then, for any $1\leq j\leq q$, we have

$\displaystyle\mathbb{E}[\widehat{\theta}_{b}^{(j)}]$
$\displaystyle=\theta+O(b^{-1})+O(T^{-1})$
$\displaystyle\mbox{Var}[\widehat{\theta}_{b}^{(j)}]$
$\displaystyle=b^{-1}\,\langle
g\,g^{\star}\,f^{2}\rangle+T^{-1}\,\langle\langle
g\,g\,F\rangle\rangle+O(T^{-2})+O(b^{-2})+O(b^{-1}T^{-1}).$

Also, for $1\leq i,j\leq q$ and $i\neq j$, we have
$\mbox{Cov}[\widehat{\theta}_{b}^{(i)},\widehat{\theta}_{b}^{(j)}]=O(T^{-1}).$

#### Proof of Theorem 2.

Let $d_{T}(\lambda)=\sum_{t=1}^{T}X_{t}e^{-it\lambda}$, so that for
$\lambda\neq 0$ we have $I_{T}(\lambda)=T^{-1}{|d_{T}(\lambda)|}^{2}$. When
evaluated at a non-zero Fourier frequency $\lambda_{\ell}$,
$d_{T}(\lambda_{\ell})$ gives the same value when computed from $X_{t}-\mu$
instead of $X_{t}$, because $\sum_{t=1}^{T}e^{-it\lambda_{\ell}}=0$. So
without loss of generality we can assume that $\mu=0$ in our analysis.
Furthermore, because
$\mathbb{E}[d_{T}(\lambda)d_{T}(-\lambda)]=\mbox{cum}(d_{T}(\lambda),d_{T}(-\lambda))$
when $\mu=0$, we can apply Theorem 4.3.2 of Brillinger (1981) to obtain

$\displaystyle\mathbb{E}[\widehat{\theta}_{b}^{(j)}]$
$\displaystyle=\frac{1}{bT}\sum_{\ell=1}^{[b/2]}\left(g(2\pi\ell/b)+g(-2\pi\ell/b)\right)\,\mathbb{E}[d_{T}(2\pi\ell/b+2\pi
j/T)d_{T}(-2\pi\ell/b-2\pi j/T)]$
$\displaystyle=\frac{1}{bT}\sum_{\ell=1}^{[b/2]}\left(g(2\pi\ell/b)+g(-2\pi\ell/b)\right)\,\left(Tf(2\pi\ell/b+2\pi
j/T)+O(1)\right)$
$\displaystyle=O(T^{-1})+\frac{1}{b}\sum_{\ell=1}^{[b/2]}\left(g(2\pi\ell/b)+g(-2\pi\ell/b)\right)\,f(2\pi\ell/b+2\pi
j/T),$

using the boundedness of $g$. The case $k=2$ of Assumption (B) implies that
$\partial_{\lambda}f(\lambda)$ is bounded in $\lambda$, and hence by a Taylor
series expansion $f(2\pi\ell/b+2\pi j/T)=f(2\pi\ell/b)+O(b^{-1})$, since
$|j|/T\leq q/T=b^{-1}$. Thus, $\mathbb{E}[\widehat{\theta}_{b}^{(j)}]$ equals
$b^{-1}\sum_{\ell=1}^{[b/2]}\left(g(2\pi\ell/b)+g(-2\pi\ell/b)\right)\,f(2\pi\ell/b)$
plus terms that are $O(T^{-1})+O(b^{-1})$. Finally,

$b^{-1}\sum_{\ell=1}^{[b/2]}\left(g(2\pi\ell/b)+g(-2\pi\ell/b)\right)\,f(2\pi\ell/b)=b^{-1}\sum_{|\ell|=1}^{[b/2]}g(2\pi\ell/b)\,f(2\pi\ell/b)$

by the evenness of $f$, and this last expression is the Riemann sum on a mesh
of size $b^{-1}$ of $\langle gf\rangle=\theta$. This proves the first
assertion.

For the variance, we write

$\displaystyle\mbox{Var}[\widehat{\theta}_{b}^{(j)}]$
$\displaystyle=\frac{1}{b^{2}T^{2}}\sum_{\ell,k=1}^{[b/2]}\left(g(2\pi\ell/b)+g(-2\pi\ell/b)\right)\left(g(2\pi
k/b)+g(-2\pi k/b)\right)$
$\displaystyle\cdot\mbox{cum}\left(d_{T}(2\pi\ell/b+2\pi
j/T)d_{T}(-2\pi\ell/b-2\pi j/T),d_{T}(2\pi k/b+2\pi j/T)d_{T}(-2\pi k/b-2\pi
j/T)\right).$

To compute the cumulant, we use Theorems 2.3.2 and 4.3.2 of Brillinger (1981).
In the case that $\ell=k$, we find the cumulant, up to terms that are $O(1)$,
is

$T^{2}{f(2\pi\ell/b+2\pi j/T)}^{2}+TF(2\pi\ell/b+2\pi j/T,-2\pi\ell/b-2\pi
j/T,2\pi\ell/b+2\pi j/T).$

Hence we obtain the “diagonal” contribution to
$\mbox{Var}[\widehat{\theta}_{b}^{(j)}]$ is

$\displaystyle\frac{1}{b^{2}T^{2}}\sum_{\ell=1}^{[b/2]}{\left(g(2\pi\ell/b)+g(-2\pi\ell/b)\right)}^{2}$
$\displaystyle\cdot\left(T^{2}{f(2\pi\ell/b+2\pi j/T)}^{2}+TF(2\pi\ell/b+2\pi
j/T,-2\pi\ell/b-2\pi j/T,2\pi\ell/b+2\pi j/T)\right),$

and the term involving the tri-spectrum is of lower order, and therefore is
vanishing asymptotically. So the the diagonal contribution is asymptotic to

$\frac{1}{2\pi
b}\int_{0}^{\pi}{(g(\lambda)+g(-\lambda))}^{2}\,{f(\lambda)}^{2}d\lambda=\frac{1}{2\pi
b}\int_{-\pi}^{\pi}g(\lambda)(g(\lambda)+g(-\lambda))\,{f(\lambda)}^{2}d\lambda=b^{-1}\langle
g(g+g^{\sharp})f^{2}\rangle.$

Turning to the “off-diagonal” contribution to the variance, we examine the
case that $\ell\neq k$, finding that the cumulant is

$TF(2\pi\ell/b+2\pi j/T,-2\pi\ell/b-2\pi j/T,2\pi k/b+2\pi j/T)$

up to terms that are $O(1)$. We can utilize a Taylor series expansion in each
of the three arguments of the tri-spectrum, finding that up to terms
$O(b^{-1})$ we have $F(2\pi\ell/b,-2\pi\ell/b,2\pi k/b)$. Hence the off-
diagonal portion of the asymptotic variance is

$\displaystyle\frac{1}{b^{2}T}\sum_{\ell\neq
k=1}^{[b/2]}\left(g(2\pi\ell/b)+g(-2\pi\ell/b)\right)\left(g(2\pi k/b)+g(-2\pi
k/b)\right)\,F(2\pi\ell/b,-2\pi\ell/b,2\pi k/b)$
$\displaystyle=\frac{{(b/2)}^{2}}{b^{2}T\pi^{2}}\int_{0}^{\pi}\int_{0}^{\pi}(g(\lambda)+g(-\lambda))(g(\omega)+g(-\omega))\,F(\lambda,-\lambda,\omega)\,d\lambda
d\omega+O(b^{-1}T^{-1}),$

using the Riemann sum approximation for both integrals. Exploiting the
symmetries of the tri-spectrum, the above expression simplifies to
$T^{-1}\langle\langle ggF\rangle\rangle+O(b^{-1}T^{-1})$.

Finally, the covariance (for $i\neq j$) is

$\displaystyle\mbox{Cov}[\widehat{\theta}_{b}^{(i)},\widehat{\theta}_{b}^{(j)}]$
$\displaystyle=\frac{1}{b^{2}T^{2}}\sum_{\ell,k=1}^{[b/2]}\left(g(2\pi\ell/b)+g(-2\pi\ell/b)\right)\left(g(2\pi
k/b)+g(-2\pi k/b)\right)$
$\displaystyle\cdot\mbox{cum}\left(d_{T}(2\pi\ell/b+2\pi
i/T)d_{T}(-2\pi\ell/b-2\pi i/T),d_{T}(2\pi k/b+2\pi j/T)d_{T}(-2\pi k/b-2\pi
j/T)\right).$

First examining the case that $\ell=k$, the cumulant involves terms that are
$O(1)$, plus

$T^{2}f(2\pi\ell/b+2\pi i/T)f(2\pi\ell/b+2\pi j/T){\bf
1}_{\\{i=j\\}}+TF(2\pi\ell/b+2\pi i/T,-2\pi\ell/b-2\pi i/T,2\pi\ell/b+2\pi
j/T).$

The indicator arises, because in the cumulant asymptotics of Theorem 4.3.2 of
Brillinger (1981), there is a Dirac operator evaluated at the sum of Fourier
frequencies $-2\pi\ell/b-2\pi i/T+2\pi\ell/b+2\pi j/T$, which is zero if and
only if $i=j$. Since this case is excluded by assumption, we find that the
diagonal contribution to the covariance is $O(T^{-1})$. For the case that
$\ell\neq k$, we now consider a Dirac operator evaluated at the sum of Fourier
frequencies $-2\pi\ell/b-2\pi i/T+2\pi k/b+2\pi j/T$, which equals zero if and
only if $j-i=q(\ell-k)$. But since $i\neq j$, we see that
$|j-i|\in\\{1,\ldots,q-1\\}$, and so the condition is impossible. Hence the
only contribution to the cumulant is the tri-spectrum, and we obtain

$O(1)+TF(2\pi\ell/b+2\pi i/T,-2\pi\ell/b-2\pi i/T,2\pi k/b+2\pi j/T).$

It follows that the covariance is $O(T^{-1}$) asymptotically, which completes
the proof. $\quad\Box$

Denote the left-hand-side of (13) by
$S_{T}(\theta)=\sqrt{T}\,\left(\widehat{\theta}_{T}-\theta\right)$; this was
previously called a ‘root’. As a consequence of Theorem 2, the $j$th oracle
skip-sample root
$S^{(j)}_{T}(\theta)=\sqrt{b}\,(\widehat{\theta}_{b}^{(j)}-\theta)$ has
asymptotic variance

$\langle g\,g^{\star}\,f^{2}\rangle+(b/T)\,\langle\langle
g\,g\,F\rangle\rangle.$ (15)

Theorem 2 allows us to take advantage of the avenue suggested by Remark 3.1,
i.e., using the skip-sampling estimator of the asymptotic variance of eq.
(13), and then use the normal tables for inference. Note that here the rate
$a_{T}=\sqrt{T}.$ Letting $b=o(\sqrt{T})$ to satisfy (8), Corollary 1 would be
applicable as long as assumption (10) were also verified. In view of Remark
3.2, the latter is expected to hold but it is cumbersome to evaluate.

However, there is an additional issue: since $b/T\rightarrow 0$, the second
term of (15) asymptotically drops out, which is undesirable in terms of
capturing the variance given in (13). To elaborate, for roots such that
$\langle\langle g\,g\,F\rangle\rangle=0$ the asymptotic variance of
$S_{T}^{(j)}(\theta)$ is correct, but otherwise must be adjusted to account
for the non-trivial contribution from the tri-spectrum. Whenever
$\langle\langle g\,g\,F\rangle\rangle=0$, we will say the asymptotic
distribution (13) is “tri-spectrum free."

###### Corollary 2

Assume the assumptions of Theorem 2, the additional assumption (10) and
$b=o(\sqrt{T})$. Then, the skip-sampling estimator $\hat{v}_{b}$ from eq. (7)
is consistent for the asymptotic variance appearing in eq. (13) when the
latter is “tri-spectrum free."

###### Remark 4.1

Requiring that the asymptotic distribution be “tri-spectrum free" is common
with several resampling methods in the frequency domain. For example, the
original frequency domain bootstrap of Franke and Härdle (1992) fails to
capture the second term of (15) even for linear processes; see Paparoditis
(2002) for a review. In general, by the Wold decomposition we have
$X_{t}=EX_{0}+\sum_{j\geq 0}\psi_{j}\epsilon_{t-j}$, where the sequence
$\epsilon_{t}$ is mean zero, uncorrelated with variance $\sigma^{2}$, i.e., a
white noise, but not necessarily independent, identically distributed
(i.i.d.). If the process $\\{X_{t}\\}$ is linear (and causal), then
$\epsilon_{t}\sim\mbox{i.i.d.}$ as well, and the expression for the variance
greatly simplifies. In this case,

$\langle\langle g\,g\,F\rangle\rangle=(\eta-3)\,{\langle g\,f\rangle}^{2},$

where $\eta=\mathbb{E}[\epsilon_{t}^{4}]/\sigma^{4}$. This yields a classical
result: for linear time series, if the innovation kurtosis is that of a
Gaussian (i.e., $\eta=3$), or in the special case when the linear spectral
mean is zero (i.e., $\langle g\,f\rangle=0$), then the asymptotic distribution
(13) is tri-spectrum free and Corollary 2 is applicable.

If $\\{X_{t}\\}$ is linear but $\eta\neq 3$, it may still be possible to
conduct inference on spectral means via a hybrid procedure employing skip-
sampling as a component. For example, let $\widehat{\eta}$ be the estimator of
$\eta$ based on the technique of Fragkeskou and Paparoditis (2016), and let
$\widehat{f}$ be a consistent estimator of the spectral density $f$. Then, we
can estimate the asymptotic variance appearing in eq. (13) by
$\widehat{v}_{b}+(\widehat{\eta}-3)\,{\langle g\,\widehat{f}\rangle}^{2}.$
Alternative hybrid methods are also available, see e.g. Janas and Dahlhaus
(1994), Kreiss and Paparoditis (2003), or Meyer et al. (2020).

###  4.2 Ratio statistics

Consider a parameter $\theta$ that is obtained as the finite ratio of two
linear spectral means, i.e., $\theta=\ \langle p\,f\rangle/\langle
m\,f\rangle$ for some fixed bounded functions $p(\lambda)$ and $m(\lambda)$
having bounded variation on $[-\pi,\pi]$. We can estimate $\theta$ by a so-
called ratio statistic given by

$\widehat{\theta}_{T}=\frac{\sum_{\ell\in
R_{T}}p(2\pi\ell/T)\,I_{T}(2\pi\ell/T)}{\sum_{\ell\in
R_{T}}m(2\pi\ell/T)\,I_{T}(2\pi\ell/T)}.$ (16)

The prime example of a ratio statistic is the sample autocorrelation at lag
$k$, where $p(\lambda)=e^{ik\lambda}$ and $m(\lambda)=1$.

Ratio statistics have an asymptotic distribution that can be tri-spectrum free
under some conditions —such as linearity of the time series— and are thus
amenable to frequency domain resampling. In fact, Dahlhaus and Janas (1996)
showed that the original frequency domain bootstrap of Franke and Härdle
(1992) is not only consistent, but higher-order accurate for ratio statistics
from linear time series. To elaborate, Remark 4.1 implies that —provided the
process is linear— ratio statistics satisfy a simplified version of (13),
namely

$\sqrt{T}\,\left(\widehat{\theta}_{T}-\theta\right)\stackrel{{\scriptstyle{\cal
L}}}{{\Longrightarrow}}\mathcal{N}\left(0,\langle
g\,g^{\star}\,f^{2}\rangle/{\langle m\,f\rangle}^{2}\right)\ \ \mbox{as}\ \
T\to\infty;$ (17)

here, $g=p-m\theta$, and the notation $g^{\star}$ was defined right after
(13). In analogy with the previous subsection, we define the $j$th skip-sample
ratio statistic via

$\widehat{\theta}_{b}^{(j)}=\frac{b^{-1}\sum_{\ell=1}^{[b/2]}p^{\star}(2\pi\ell/b)\,I_{T}(2\pi\ell/b+2\pi
j/T)}{b^{-1}\sum_{\ell=1}^{[b/2]}m^{\star}(2\pi\ell/b)\,I_{T}(2\pi\ell/b+2\pi
j/T)}.$

###### Theorem 3

Assume that $\\{X_{t}\\}$ is a strictly stationary linear process that
satisfies Assumption B. Consider a finite ratio of linear spectral means
$\theta=\langle p\,f\rangle/\langle m\,f\rangle$, for some fixed bounded
functions $p(\lambda)$ and $m(\lambda)$ having bounded variation, and let
$\widehat{\theta}_{T}$ be the ratio statistic (16). Let $b$ be a sequence
satisfying (6), and set $g=p-m\theta$. Then, for any $1\leq j\leq q$, we have

$\displaystyle\mathbb{E}[\widehat{\theta}_{b}^{(j)}]$
$\displaystyle=\theta+O(b^{-1})+O(T^{-1})$
$\displaystyle\mbox{Var}[\widehat{\theta}_{b}^{(j)}]$
$\displaystyle=b^{-1}\,\langle g\,g^{\star}\,f^{2}\rangle/{\langle
m\,f\rangle}^{2}+O(T^{-2})+O(b^{-2})+O(b^{-1}T^{-1}).$

Also, for $i\neq j$ and $1\leq i,j\leq q$, we have
$\mbox{Cov}[\widehat{\theta}_{b}^{(i)},\widehat{\theta}_{b}^{(j)}]=O(T^{-1}).$

#### Proof of Theorem 3.

It is easy to show that

$\widehat{\theta}_{b}^{(j)}-\theta=\frac{b^{-1}\sum_{\ell=1}^{[b/2]}g^{\star}(2\pi\ell/b)\,I_{T}(2\pi\ell/b+2\pi
j/T)}{b^{-1}\sum_{\ell=1}^{[b/2]}m^{\star}(2\pi\ell/b)\,I_{T}(2\pi\ell/b+2\pi
j/T)}.$

The denominator, denoted by $\widehat{\langle m\,f\rangle}$ for short,
converges in probability to $\langle m\,f\rangle$ by Theorem 2. Because
$\theta$ is finite, this limit must be non-zero. Letting the numerator be
denoted $\widehat{\langle g\,f\rangle}$, we obtain

$\widehat{\theta}_{b}^{(j)}-\theta=\frac{\widehat{\langle
g\,f\rangle}}{\langle m\,f\rangle}+\frac{\widehat{\langle
g\,f\rangle}\left(\langle m\,f\rangle-\widehat{\langle
m\,f\rangle}\right)}{\langle m\,f\rangle\,\widehat{\langle m\,f\rangle}},$

and the second term has mean that is $O(b^{-1})+O(T^{-1})$, using the Cauchy-
Schwarz inequality, the delta method, and the variance results of Theorem 2.
This is because the mean of $\widehat{\langle g\,f\rangle}$ is zero plus lower
order terms. Finally, we can compute the mean and variance of
$\widehat{\langle g\,f\rangle}/\langle m\,f\rangle$ using Theorem 2, and
obtain the stated results, noting that $\langle\langle ggF\rangle\rangle=0$
since the process is linear and $\langle g\,f\rangle=0$. $\Box$

Theorem 3 confirms the validity of assumption set (9) in the context of ratio
statistics; the following corollary then ensues.

###### Corollary 3

Assume the assumptions of Theorem 3, the additional assumption (10) and
$b=o(\sqrt{T})$. Then, the skip-sampling estimator $\hat{v}_{b}$ from eq. (7)
is consistent for the asymptotic variance appearing in eq. (17).

## Acknowledgments

This report is released to inform interested parties of research and to
encourage discussion. The views expressed on statistical issues are those of
the authors and not those of the U.S. Census Bureau. Research of the second
author partially supported by NSF grant DMS 19-14556.

## References

  * [1] Brillinger, D.R. (1981). Time Series: Data Analysis and Theory, Holden-Day, New York.
  * [2] Bertail, P. and Dudek, A.E. (2021). Consistency of the frequency domain bootstrap for differentiable functionals, Electronic Journal of Statistics 15, 1–36.
  * [3] Brockwell, P. J. and Davis, R. A. (1991). Time Series: Theory and Methods, 2nd ed., Springer, New York.
  * [4] Dahlhaus, R. (1985). Asymptotic normality of spectral estimates, Journal of Multivariate Analysis 16, 412–431.
  * [5] Dahlhaus, R. and Janas, D. (1996). A frequency domain bootstrap for ratio statistics in time series analysis, Ann. Statist. 24(5), 1934–1963.
  * [6] Efron, B. (1979). Bootstrap methods: another look at the jackknife, Ann. Statist. 7, 1–26.
  * [7] Franke, J. and Härdle, W. (1992). On bootstrapping kernel spectral estimates, Ann. Statist. 20, 121–145.
  * [8] Fragkeskou, M. and Paparoditis, E. (2016). Inference for the fourth-order innovation cumulant in linear time series, Journal of Time Series Analysis 2(37), 240–266.
  * [9] Hurvich, C. M. and Zeger, S. L. (1987). Frequency domain bootstrap methods for time series, New York University Working Paper.
  * [10] Janas, D. and R. Dahlhaus (1994). A frequency domain bootstrap for time series. In: Computationally Intensive Statistical Methods. Proceedings of the 26th Symposium on the Interface, (J. Sall and A. Lehman, eds.). Interface Foundation of North America, Faifax Station, VA, pp. 423–425
  * [11] Kirch, C. and D.N. Politis (2011). TFT-bootstrap: resampling time series in the frequency domain to obtain replicates in the time domain, Ann. Statist. 39(3), 1427–1470.
  * [12] Kreiss, J.-P. and Paparoditis, E. (2003). Autoregressive aided periodogram bootstrap for time series, Ann. Statist. 31(6), 1923–1955.
  * [13] Kreiss, J.-P. and Paparoditis, E. (2023). Bootstrap for Time Series: Theory and Applications, Springer, Heidelberg.
  * [14] Lahiri, S.N. (2003a). Resampling Methods for Dependent Data, Springer, New York.
  * [15] Lahiri, S.N (2003b). A necessary and sufficient condition for asymptotic independence of discrete Fourier transforms under short-and long-range dependence, Ann. Statist. 31(2), 613–641.
  * [16] McElroy, T.S. and Politis, D.N. (2020). Time Series: A First Course with Bootstrap Starter, Chapman and Hall/CRC Press, Boca Raton.
  * [17] McElroy, T.S. and Roy, A. (2022). Model identification via total Frobenius norm of multivariate spectra, Journal of the Royal Statistical Society, Series B 84, 473–495.
  * [18] Meyer, M. Paparoditis, E. and Kreiss, J.-P. (2020). Extending the validity of frequency domain bootstrap methods to general stationary processes, Ann. Statist. 48(4), 2404–2427.
  * [19] Paparoditis, E. (2002). Frequency domain bootstrap for time series. In Empirical Process Techniques for Dependent Data (H. Dehling et al., eds.) 365–381. Birkhäuser, Boston.
  * [20] Paparoditis, E. and Politis, D.N. (1999). The local bootstrap for periodogram statistics, Journal of Time Series Analysis 20, 193–222.
  * [21] Politis, D.N., Romano, J.P. and Wolf, M. (1999). Subsampling, Springer, New York.
  * [22] Taniguchi, M. and Kakizawa, Y. (2000). Asymptotic Theory of Statistical Inference for Time Series, Springer, New York.


# Developing an Automated Detection, Tracking and Analysis Method for Solar
Filaments Observed by CHASE via Machine Learning

Z. Zheng School of Astronomy and Space Science, Nanjing University, Nanjing
210023, China  Key Laboratory of Modern Astronomy and Astrophysics, Ministry
of Education, Nanjing 210023, China  Q. Hao School of Astronomy and Space
Science, Nanjing University, Nanjing 210023, China  Key Laboratory of Modern
Astronomy and Astrophysics, Ministry of Education, Nanjing 210023, China  Y.
Qiu Institute of Science and Technology for Deep Space Exploration, Suzhou
Campus, Nanjing University, Suzhou 215163, China  J. Hong School of Astronomy
and Space Science, Nanjing University, Nanjing 210023, China  Key Laboratory
of Modern Astronomy and Astrophysics, Ministry of Education, Nanjing 210023,
China  C. Li School of Astronomy and Space Science, Nanjing University,
Nanjing 210023, China  Key Laboratory of Modern Astronomy and Astrophysics,
Ministry of Education, Nanjing 210023, China  Institute of Science and
Technology for Deep Space Exploration, Suzhou Campus, Nanjing University,
Suzhou 215163, China  M.D. Ding School of Astronomy and Space Science, Nanjing
University, Nanjing 210023, China  Key Laboratory of Modern Astronomy and
Astrophysics, Ministry of Education, Nanjing 210023, China Q. Hao, C. Li
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Studies on the dynamics of solar filaments have significant implications for
understanding their formation, evolution, and eruption, which are of great
importance for space weather warning and forecasting. The H$\alpha$ Imaging
Spectrograph (HIS) onboard the recently launched Chinese H$\alpha$ Solar
Explorer (CHASE) can provide full-disk solar H$\alpha$ spectroscopic
observations, which bring us an opportunity to systematically explore and
analyze the plasma dynamics of filaments. The dramatically increased
observation data require automate processing and analysis which are impossible
if dealt with manually. In this paper, we utilize the U-Net model to identify
filaments and implement the Channel and Spatial Reliability Tracking (CSRT)
algorithm for automated filament tracking. In addition, we use the cloud model
to invert the line-of-sight velocity of filaments and employ the graph theory
algorithm to extract the filament spine, which can advance our understanding
of the dynamics of filaments. The favorable test performance confirms the
validity of our method, which will be implemented in the following statistical
analyses of filament features and dynamics of CHASE/HIS observations.

Solar filaments — Convolutional neural networks — Astronomy image processing

##  1 Introduction

Solar filament is one of the typical solar activities in solar atmosphere,
which is about 100 times cooler and denser than its surrounding corona
(Labrosse et al., 2010). They are observed as dark elongated structures with
several barbs, but are seen as bright structures suspended over the solar limb
called prominences (Vial & Engvold, 2015). Filaments are always align with
photospheric magnetic polarity inversion line (PIL), where the magnetic flux
cancellation often takes place (Martin, 1998; Vial & Engvold, 2015). Filaments
sometimes undergo large-scale instabilities, which break their equilibria and
lead to eruptions. There is a close relationship among the erupting filaments,
flares, and coronal mass ejections, which are different manifestations of one
physical process at different evolutionary stages (Gopalswamy et al., 2003).
Therefore, the study of the formation, evolution and eruption of filament is
not only of great significance to understand the essence physics of solar
activities, but also of practical significance for accurately predicting the
hazardous space weather (Chen, 2011; Chen et al., 2020).

Filament are usually observed by ground-based solar H$\alpha$ telescopes
around the world, such as Meudon, Big Bear, Kanzelhöhe, Kodaikanal, and
Huairou. These telescopes have been the work horses of most of the current
knowledge on filaments (Chatzistergos, Theodosios et al., 2023). To study the
mechanisms of solar eruptions and the plasma dynamics in the lower atmosphere,
the Chinese H$\alpha$ Solar Explorer (CHASE; Li et al., 2019, 2022) was
launched into a Sun-synchronous orbit on October 14, 2021. The scientific
payload onboard CHASE is the H$\alpha$ Imaging Spectrograph (HIS; Liu et al.,
2022), which can provide solar H$\alpha$ spectroscopic observations. It brings
us an opportunity to systematically explore and analyze the plasma dynamics of
filaments in details. At the same time, the data volume of CHASE/HIS
observations has dramatically increased, which also brings challenges to
efficiently process such huge amount of data.

In order to statistically obtain the filament features, Gao et al. (2002)
developed an automated algorithm combining the intensity threshold and the
region growing method. Since then, a number of automated filament detection
methods and algorithms based on classical image processing techniques have
been developed in the past decades (Shih & Kowalski, 2003; Fuller et al.,
2005; Bernasconi et al., 2005; Qu et al., 2005; Wang et al., 2010; Labrosse et
al., 2010; Yuan et al., 2011; Hao et al., 2013, 2015). Shih & Kowalski (2003)
adopted local thresholds which were chosen by median values of the image
intensity to extract filaments. However, this kind of threshold selection
cannot guarantee robust results since the bright features on images can
significantly affect the value of the thresholds. To overcome this problem,
some authors have developed the adaptive threshold methods (Qu et al., 2005;
Yuan et al., 2011; Hao et al., 2015). Particularly, Qu et al. (2005) applied
the Support Vector Machine (SVM) technique to distinguish filaments from
sunspots. The development of graphics processing unit (GPU) and machine
learning in recent years brings a powerful set of techniques which drive
innovations in areas of computer vision, natural language processing,
healthcare, and also in astronomy in the recent decade (Smith & Geach, 2023;
Asensio Ramos et al., 2023). Artificial neural networks, especially the
convolutional neural networks (CNNs) have been leading the trend of machine
learning on the feature segmentation for years since AlexNet (Krizhevsky et
al., 2012). Recently, CNNs are widely adopt to automated detect filament, and
they have been proven to have a high performance (Zhu et al., 2019; Liu et
al., 2021; Guo et al., 2022). Zhu et al. (2019) developed an automated
filament detection method based on the U-Net (Ronneberger et al., 2015), a
kind of deep learning architecture,which has an excellent performance in
semantic segmentation since it can gain strong robustness even from a small
data set. Their test performance showed that the new method can segment
filaments directly and avoid generating segmentation with a large number of
noise points as classical image processing method. A filament may be split
into several fragments during its evolution. To figure out whether they belong
to one filament is crucial for the study of filament evolution. Many authors
adopted morphological operations, the distance criterion and the slopes of the
fragments (Shih & Kowalski, 2003; Fuller et al., 2005; Bernasconi et al.,
2005; Qu et al., 2005; Hao et al., 2013, 2015). However, these thresholds for
distance or angles do not always work. Guo et al. (2022) proposed a new method
base on deep-learning instance-segmentation model CondInst (Tian et al., 2020,
2023) to solve the fragment problems. Since such methods are supervised
machine learning methods, we have to first provide data set with filament
labeled. The labeled data are usually obtained by manually annotating the
ground-based H$\alpha$ images, which cannot maintain the consistence and
accuracy of labelling. Thanks to the CHASE mission for providing the seeing-
free H$\alpha$ spectroscopic observations, we can get the precise boundaries
of filaments for training the deep learning models.

In this paper, we developed an efficient and robust automated detection and
tracking method for filament observed by CHASE. Figure 1 shows the flowchart
of our method consisting of three parts: data preparation, filament detection,
and filament tracking. In Section 2, we describe the data preparation,
including the filament labeling, calibration for image data, and other
necessary adjustments. The pipelines for the automated detection and tracking
system are descried in Section 3 and 5, followed by the description of
performance, respectively. Section 4 is dedicated to the inversion of line-of-
sight velocity and filament spine extraction. Discussion and conclusions are
given in Section 6.

Figure 1: The flowchart of our data preparation, detection, and tracking
methods. Solid arrows represent the flow of data processing, while dashed
arrows denote data connections between different stages of the processing.

##  2 Data preparation

Figure 2: An observation obtained on 2023 January 19 as an example to show the
schematic of data preparation. (a) The original H$\alpha$ line center image.
(b) The image after limb-darkening removal. (c) The K-Means spectral
classification of filament candidates with yellow contours. (d) An enlarged
part of (c). (e) and (f) is similar to (c) and (d), but is the final labeled
filament after area filtering and morphological close operation.

The CHASE/HIS implements raster scanning of the full solar disk within 60
seconds at the wavebands of H$\alpha$ and Fe I with a spectral sampling of
0.048 Å and a pixel resolution of $1.04^{\prime\prime}$ in the binning mode
(Qiu et al., 2022). Although our detection and tracking method are based on
the H$\alpha$ line center images, the high quality H$\alpha$ spectral
information allows us to get the precise boundaries of filaments by spectral
classification for training the deep learning models. As shown in Figure 1,
before the filament detection and tracking, we proposed a series of
preprocessings to build a data set for training and testing the deep-learning
model. First, we used a second-order cosine function to remove the limb
darkening. Next, we used a unsupervised learning method called K-means
(MacQueen, 1967) to classify the whole disk spectra in order to get the
precise boundaries of filaments. The following subsections explain the details
of each approach.

###  2.1 Spectral Classification for Filament Labeling

Our filament detection method is a supervised machine learning method, which
means we need to provide the data sets with filament labeled. Here we utilize
the CHASE full-disk H$\alpha$ spectra to identify and label the filament. It
is difficult to identify filaments by spectra manually since there are 118
wavelength points in CHASE H$\alpha$ profile. We employed K-Means algorithm
(MacQueen, 1967), an unsupervised clustering algorithm which has been widely
applied in spectral classification of solar physics (Viticchié & Sánchez
Almeida, 2011; Panos et al., 2018; Asensio Ramos et al., 2023). It requires
the initial setting of the number of clusters, denoted as $k$. Subsequently,
it will assign $k$ centroids $(\mu_{1},\mu_{2},\cdots,\mu_{k})$ and divide the
data into $k$ clusters $(S_{1},S_{2},\cdots,S_{k})$ based on the distance to
the centroids. Then, it continually updates the centroids and clusters to
minimize the intracluster distance by solving following equation:

$\mathop{\arg\min}_{\mu}\sum_{i=1}^{k}\sum_{x_{j}\in S_{i}}\|x_{j}-\mu_{i}\|,$
(1)

where $x$ is the contrast of the spectral intensities of each pixel with the
average spectral intensity $I/I_{avg}$.

The K-means method is very sensitive to the uneven radiation intensity
distribution across the solar disk, which is adversely affecting spectral
classification for labeling filaments. We use a second-order cosine function
(Pierce & Slaughter, 1977) to remove the limb darkening:

$I_{\lambda}^{*}(R)=a_{\lambda}+b_{\lambda}*cos\theta+c_{\lambda}*cos^{2}\theta,$
(2)

where $I_{\lambda}^{*}(R)$ is the mean observed radiation intensity at the
radius $R$ in wavelength $\lambda$, $R_{S}$ is the radius of the Sun,
$cos\theta=\sqrt{1-(R/R_{s})^{2}}$ and $a_{\lambda}$, $b_{\lambda}$,
$c_{\lambda}$ are the parameters to be fitted. We used the radiation intensity
along the solar equator as the fitting data and applied a least-squares fit to
minimize the deviation rate. The radiation intensity along the solar equator
is selected since there are few activities and it is longer enough to fit the
whole disk. For each wavelength $\lambda$ of the H$\alpha$ profile, we need
find the best $a_{\lambda}$, $b_{\lambda}$, $c_{\lambda}$ to minimize the
following equation:

$\sum_{p_{i}\in eq}\mid
I_{\lambda}(p_{i})-I_{\lambda}^{*}(R(p_{i}))\mid/I_{\lambda}^{*}(R(p_{i})),$
(3)

where $I_{\lambda}^{*}(p_{i})$ is the observed radiation intensity at point
$p_{i}$ of the equator and $R(p_{i})$ is the radius at $p_{i}$. We limited the
fitting to 880 pixels (about $0.95R_{s}$) since the pixels alone the edge may
vary with the observation. Figure 2(a) and (b) give an example to show the
H$\alpha$ line center image before and after limb-darkening removal,
respectively.

After a series of trial and error we set the $k=30$, i.e., the K-Means
algorithm automated categorizes the spectra into 30 classes, in which they
belong to sunspot, plage, filament and so on. Then we manually select the
class that most closely matches the filament as the output label candidate. As
shown in Figure 2(c) and (d), the K-Means spectral classification could
effectively segment filaments from the solar disk. However, small
chromospheric fibers can not be removed since they are also cold material and
have similar H$\alpha$ spectral profiles to the filaments. Here, we set a area
threshold of 64 pixels ( about 69 arcsec2) to sieve the chromospheric fibers.
We can also find that the filament classified by K-Means are more mottled with
holes as shown in Figure 2(d). We adopt the morphological close operation with
round structure element (radius is 5 pixel) to fill these holes. Figure 2(e)
and (f) show the results after adjustment. Then the labeled data are used for
model training. In addition, we found that the data with morphological close
operation can improve the ability of U-Net model to distinguish the sunspot
and filament. It is because the input of U-Net model is the normalized
H$\alpha$ line center images, which can not distinguish the sunspot and
filament only by intensity since they have similar dark appearance.

###  2.2 Data arrangement for model training

We collected 120 sets of H$\alpha$ spectral observation from CHASE and
constructed our labeled dataset by applying the approaches mentioned above.
The data span from December 2022 to July 2023, with a time interval of about 2
days. In a supervised machine learning method, the data usually are divided
into training, validation, and testing set, respectively. The validation set
was utilized during the training phase to choose models, which can prevent
overfitting on the training set. The testing set was employed to assess the
performance of trained models. Therefore, these 120 sets are divided into 20
groups, the third and the sixth sets in each group are chosen as the
validation and testing sets, respectively, while the others are chosen as the
training set, with the ratio being about $1:1:4$. In this way, we can ensured
a evenly distribution in our data sets so that the trained model could be
robust and increase its generalization performance.

The data augmentation technique is usually applied to enlarge the dataset by
flipping and rotating transformations since the convolution kernel does not
have rotational symmetry and axis-symmetry, which can reduce the risk of
overfitting (Mikołajczyk & Grochowski, 2018). Here, we also adopt data
augmentation during the training process to enhance the robustness of our
model. In each training iteration, the data are randomly applied
transformations, including vertical flipping (with a $50\%$ probability) and
rotation with a angle within $[-45^{\circ},45^{\circ}]$ (with an $80\%$
probability).

Figure 3: Schematic representation of our U-Net model. The model takes
H$\alpha$ line-center images as input and produces a binary image of the same
size as output. The cubes represent feature maps, where the dimension of each
map is indicated on its left, and the number of channels is indicated above
it. Operations for each channel are represented by arrows with different
colors.

##  3 Filament detection

Recently, CNNs are adopt to automated detect filament, and it have been proven
to have a high performance (Zhu et al., 2019; Liu et al., 2021; Guo et al.,
2022; Zhang et al., 2023). Here we employ the U-Net (Ronneberger et al., 2015)
architecture, a kind of CNN architecture, to implement filament detection in
our work. Figure 1 gives the flowchart of the detection method. Compared to
using H$\alpha$ spectral data as the input of K-means method during the
preprocessing steps, we only adopt the H$\alpha$ line center image as the
input of U-Net model. Note that the input images are normalized by dividing
its mean intensity instead of the maximum intensity since the errors arising
from variations in observation times, especially during the flaring time.

###  3.1 U-Net architecture for filament detection

The U-Net model derives its name from its U-shaped architecture, as shown in
Figure 3. Upon receiving image inputs, it initiates an encoding process to
progressively extract high-level semantic information while reducing the image
dimensions, as shown in left part of Figure 3. Following this, the decoding
process reconstructs the feature maps to the original image size (right part
of the U-shaped architecture), delivering detailed pixel-to-pixel level
prediction results.

In our model, the encoding and decoding parts each has four blocks, which are
repeating patterns of layers. During the encoding process, each block has two
$3\times 3$ convolution layers, followed by a $2\times 2$ max-pooling layer,
as indicated by the blue, red, and black arrows in the left side of Figure 3.
A convolution layer works as a feature detector to get the feature maps. As
the layers of the neural network go deeper, the number of the output channels
gets larger, which can regard as the feature numbers increasing. Max-pooling
layers output a feature map containing the most prominent features of the
previous feature map, which play a critical role in compressing the size of
feature maps while preserving important features and relationships in the
input images. For example, the input is an image with 1 channel and dimension
size of $2048\times 2048$, the output after the first block becomes a feature
map with 64 channels and dimension size of $1024\times 1024$, as shown in
Figure 3. The feature maps from the convolution layer passes on to activation
functions to make the network adapt to nonlinear features from the previous
layers. We use LeakyReLU and Sigmoid functions as our activation functions,
which are defined as:

$LeakyReLU_{\alpha}(x)=\left\\{\begin{array}[]{cc}\alpha x,&{x<0}\\\ x,&{x\geq
0}\end{array};\right.$ (4)

$Sigmoid(x)=1/(1+e^{-x}),$ (5)

where $\alpha$ is a hyperparameter and set to $0.01$ in our model. The Sigmoid
function is only applied in the output layer of the model to transform the
output into the range between 0 and 1 for binary classification, enabling the
distinction between targets, i.e., filaments and the background.

In the decoding process, the feature and spatial information through a series
of transpose convolutions return to the original image dimension. Transpose
convolution, often referred to as deconvolution or up-sampling, involves the
introduction of zero values through interpolation in the image, followed by a
convolution operation, allowing for effective image magnification. After this
block, the output feature map size are 4 times larger, from $128\times 128$ to
$256\times 256$, $256\times 256$ to $512\times 512$ and so on, as indicated by
Figure 3. Copy and crop, or skip connections, indicated by the gray arrows in
Figure 3, entails the concatenation of semantic information derived from the
encoding process with the corresponding feature maps during the decoding
process. This operation facilitates the transmission of augmented semantic
information, thereby bolstering the segmentation performance of the model.

Loss functions play an important role in machine learning. They define an
objective by which the performance of the model is evaluated. The parameters
learned by the model are determined by minimizing a chosen loss function. In
other words, loss functions are a measurement of how good our model is at
segmenting the filaments. We utilize the Focal Loss function (Lin et al.,
2017), a well-established choice in binary classification scenarios:

$\displaystyle FocalLoss_{\gamma,w}(y,\hat{y})=\sum-
wy(1-\hat{y})^{\gamma}\ln\hat{y}$
$\displaystyle-(1-y)\hat{y}^{\gamma}\ln(1-\hat{y}),$ (6)

where $y$ represents the values in the label, and $\hat{y}$ is the probability
values output by the model. $\gamma$ is a hyperparameter for reducing the
relative loss for well-classified examples, putting more focus on hard,
misclassified examples, i.e., filament regions. $w$ stands for the weight of
filament regions. Focal Loss proves beneficial in addressing the challenges
posed by imbalanced distribution of positive and negative samples,
corresponding to the scenario where the filament regions are much smaller than
the non-filament regions. The training of the U-Net model involves the
iterative adjustment of the weights of various convolution kernels within the
convolution layers to minimize the Focal Loss by the optimizer, which we
choose ADAM optimizer (Kingma & Ba, 2014) in our model. A dropout layer is a
regularization technique employed during model training to stochastically zero
out a subset of weights, which can mitigate overfitting in the model. We also
adapt dropout layers in our model. The detailed parameter settings are list in
Table 1.

Table 1: Parameters of Our U-Net Model. Convolution Kernels | See Figure 3
---|---
Active functions | See Figure 3
Optimizer | ADAM
Epochs | 500
Batchsize | 1
Learning rate | 1.0$\times 10^{-5}$
$\alpha$ in LeakyReLU | 0.1
Dropout rate | 0.2
$w$ in Focal Loss | 2
$\gamma$ in Focal Loss | 4

###  3.2 Performance of detection

The precision $P$ and recall ratio $R$ are common semantic segmentation
evaluation strategy which we also adopt in our model evaluation. These two
metrics are defined as:

$\displaystyle P=\frac{TP}{TP+FP}\,,$ (7) $\displaystyle
R=\frac{TP}{TP+FN}\,,$ (8)

where TP, FP, and FN denote the true positive, false positive, and false
negative measurements, respectively. In addition, the intersection over union
(IoU) is often used to evaluate the performance of a model (Rezatofighi et
al., 2019), which is defined as:

$\text{IoU}=\frac{TP}{TP+FP+FN}=\frac{R*P}{R+P-R*P}\,\,.$ (9)

In practical, the IoU over $60\%$ is a good enough score. The average
precision and recall ratios of the training, validation, and testing are
listed in Table 2. The precision, recall, and IoU ratios of each data file are
also plotted in Figure 5. These results indicated our model is a viable
strategy for solar filament recognition.

Table 2: Performance of Our U-Net Model. | IoU | Precision | Recall
---|---|---|---
Training set | 0.69 | 0.79 | 0.85
Validation set | 0.67 | 0.78 | 0.84
Testing set | 0.67 | 0.79 | 0.83
Figure 4:  The detected results by our U-Net model. (a) The detected filaments
with blue contours are plotted above the original H$\alpha$ line center
images. (b) A enlarged part of (a), where the ground-truth indicated by the
yellow contours.  Figure 5:  The performance of our U-Net model. (a), (b), and
(c) show the precision, recall, and IoU ratios of the samples in the training,
validation, and testing set, respectively. Figure 6:  A detected filament of
the CHASE observation obtained on 2023 April 25 as an example to show the
schematic of filament feature extraction. (a) The detected filament with blue
contour, between the red and blue contours are the region set as the quiet
region for LOS velocity inversion. (b) The LOS velocity distribution of the
filament. (c) The filament skeleton marked with yellow curve. (d) The filament
spine marked with orange curve. (e) The straightened filament. (f) The
corresponding LOS velocity distribution of (e). (g) The scatter plot of LOS
velocity perpendicular to the filament spine.

##  4 Filament Feature Extraction

After the filament detection, we can obtain the basic morphological features
such as filament location and area. These morphological features are not
enough for the analysis of the evolution and eruption of filaments, so we
apply the cloud model to invert line-of-sight (LOS) velocities of filaments
and employed the graph theory algorithm to extract the filament spines, which
are explained in following two subsections.

###  4.1 LOS velocity inversion of filaments

The dynamic evolution of filaments has consistently been a key issue in the
study of filament (Chen et al., 2020). CHASE provides the full-disk H$\alpha$
spectral data, which bring us a chance to extract the dynamic information of
filaments. We applied the cloud model (Beckers, 1964) to fit the H$\alpha$
spectra of filaments in order to derive its LOS velocities. The contrast
function of cloud model refers to:

$C(\lambda)=\frac{I_{0}(\lambda)-I(\lambda)}{I_{0}(\lambda)}=(1-\frac{S}{I_{0}(\lambda)})(1-e^{(-\tau)}),$
(10)

where $\tau=\tau_{0}e^{-(\lambda-\lambda_{0}-v\lambda_{0}/c)^{2}/w^{2}}$, $I$
represents the spectral profile of the filament region, $I_{0}$ corresponds to
the quiet region, $\lambda_{0}$ is the line center wavelength, and $c$ is the
light speed. The parameters to be fitted include the LOS velocity of the cloud
$v$, optical depth $\tau_{0}$, source function $S$, and Doppler width $w$.

The cloud model inversion is applied to derive the spectral profiles of each
individual filament. As shown in Figure 6(a), the blue contour represents the
filament detected by our U-Net model. The red contour represents the extension
of the blue one by 10 pixels. The average spectrum between the red and green
contours is extracted as $I_{0}$. we determine the line center wavelength
$\lambda_{0}$ by using the moment method (Yu et al., 2020). Then, the contrast
profile of the spectrum in blue contour is applied by a least-squares fitting,
providing the inversion result of the spectral profile. Figure 6(b) presents
the inversion result of LOS velocity distribution of the entire filament
region.

###  4.2 Filament spine extraction

Filament spine defines the skeleton of a filament along its full length, with
several extend barbs. In order to derive the filament spine, many authors
employed iterations of the morphological thinning and spur removal operations
(Fuller et al., 2005; Qu et al., 2005; Hao et al., 2013). However, after
iterated spur removal operations the spines often become shorter than the
original ones. Bernasconi et al. (2005) used the Euclidean distance method to
overcome the shortcoming by morphological spur operation. Yuan et al. (2011)
and Hao et al. (2015, 2016) used the algorithm based on the graph theory by
calculating the paths from end to end points of filament skeleton in pixels,
where the longest path is kept as the spine. This method is also employed in
our work.

Following the acquisition of the LOS velocity field distribution of filaments,
we want to systematically analyze the evaluations of various filaments.
However, different filaments have different length, we need compare them in a
standardized manner. Thus, we straighten and normalize the filaments along
their spines. This approach can help to uncover unified patterns in the
dynamical evolution of various kind of filaments. After extracting the spine
of the filament, we could obtain the distribution of the coordinates
$\\{(x_{i},y_{i})\\}_{i=1}^{n}$ of points along the spine with respect to the
distance set $\\{l_{i}\\}_{i=1}^{n}$, where $n$ is the number of pixels along
the spine. Subsequently, we perform spline interpolation then obtain the
parametric equation of the main spine:

$\left\\{\begin{array}[]{cc}x=X(l),\\\ y=Y(l).\end{array}\right.$ (11)

Then, we differentiate the parametric equations and apply Gaussian filter
(with $\sigma=5$) for smoothing. With the tangent equation
$(x^{\prime},y^{\prime})=(\mathrm{d}X/\mathrm{d}l,\mathrm{d}Y/\mathrm{d}l)$,
we can obtain the corresponding normal equation $(y^{\prime},-x^{\prime})$.
Subsequently, we shift each point on the major axis in the direction
perpendicular to the normal by the size of $2S_{f}/L_{f}$, where $S_{f}$ is
the filament area and $L_{f}$ is the spine length, respectively. This process
can straighten filaments along their spines according to their irregular
shapes. Figure 6 gives an example of our process. The yellow and orange curves
shows the derived filament skeleton and spine in Figure 6(c) and (d),
respectively. Figure 6(e) shows the straightened filament, which also
effectively preserves the morphological characteristics of filaments, such as
brab structures.

Moreover, the straightening approach enables a quantitative study of the
distribution of physical information along and perpendicular to the spine. As
shown in Figure 6(f) and (g), the LOS velocity distribution along and
perpendicular to its spine. Note that the majority of filaments exhibit a
rhombic or elliptical LOS velocity distribution, suggesting that the locations
with the maximum LOS velocity in filaments are typically near the spine.

Figure 7: An example for explaining our filament tracking method. (a) The
previous frame H$\alpha$ line center images observed at 2023 September 16
04:16 UT. The orange, blue, and yellow boxes with ID numbers indicate the
tracking results in different steps. (b) Similar to (a), but without the
background image. The black regions are the detected filaments by our U-Net
model. (c) The tracking results of the subsequent frame observed at 2023
September 16 05:51 UT by the CSRT tracker. (d) The results of the subsequent
frame after the initialized step. Each red box represent a merged filament.
(e) The final results after update step of the subsequent frame. (f) Similar
to (e), but with H$\alpha$ line center images as its background. Figure 8:  An
example showing the tracking results. (a) The first frame of H$\alpha$ line
center image. The orange boxes with numbers indicate the tracking results.
(b–p) A series of filament tracking results. Each panel has the same field-of-
view as the blue box in (a), and the top sub-panel is the H$\alpha$ line
center image and the bottom one is the corresponding LOS velocity
distribution.

##  5 Filament Tracking

Automated filament tracking assumes a pivotal role in statistical analysis of
filament evolution. Hao et al. (2013) considered the detected filament
locations, areas, and the differential rotations to trace the daily filament
evolution. Bonnin et al. (2013) first retrieved the location and morphology of
filament main skeleton by the filament automated recognition method by Fuller
et al. (2005), then applied a curve-matching algorithm to determine the
filaments in different frames whether they are the same one. Actually,
filament sometimes split into several fragments or partially erupted, which
represent intricate and dynamic evolution (Shen et al., 2012; Liu et al.,
2012; Sun et al., 2023; Hou et al., 2023). Tracking filament by extracting and
parameterizing certain morphological features, is only effective for
relatively larger filament and shorter time intervals. Here, we propose the
Channel and Spatial Reliability Tracking (CSRT) algorithm without requiring
additional feature extraction of transformation manually. The CSRT method is
proficient in tracking moving and deformable targets, which is quite suitable
for filament tracking.

###  5.1 Tracking method

The CSRT tracker is C++ implementation of the CSR-DCF (Channel and Spatial
Reliability of Discriminative Correlation Filter) tracking algorithm (Lukežic
et al., 2017) in OpenCV library (Bradski, 2000). The tracked object is
localized by summing the responses of the learned correlation filters and
weighted by the estimated channel reliability scores. In other words, the CSRT
tracker distinguish the target and background based on adjusting the weights
of different channels, e.g., 10 HoGs (histogram of oriented gradients)
channels and intensity channel for grayscale-images. After that, the target is
localized by the probability map and its region is updated. Furthermore,
Farhodov et al. (2019) integrated the region-based CNN pre-trained object
detection model and the CSRT tracker, and got better tracking results since
the detection model has already separated the traget and the background.
Therefore, we integrate the U-Net model and the CSRT algorithm, i.e., the
filaments detected by the U-Net model are used as inputs of the CSRT tracker.

The third part of Figure 1 shows our tracking processing scheme, which
consists of two part, i.e., the initialization step and the update step.
During the initialization step, we need to set a series of boxes which contain
the detected filaments in each frame since the input of the CSRT tracker is an
image with targets with bounding boxes. A simple way is using the minimum
boundary boxes of the detected filaments directly as the input. However,
filament sometimes may split as several fragments, or several fragments may
merge into one filament. We adopt a distance criterion of 50 pixels (about
52′′) to combine filament fragments and a area criterion of 200 pixels (about
216 arcsec2) to filter relative small filament fragments. These operation can
enhance the stability and accuracy of the CSRT tracker. Figure 7 shows two
frames of H$\alpha$ line center images observed at 04:16 UT and 05:51 UT on
2023 September 16 as the previous and subsequent frames, respectively. The
orange boxes in Figure 7(a) and (b) are the results after the initialization
step, where each box indicates the merged single filament. The previous frame
has been marked with a unique tracking ID, as shown by the numbers above the
boxes in Figure 7(a) and (b). Then the CSRT tracker tracks the subsequent
frame and output the tracked filaments based on the previous tracking IDs.
Figure 7(c) shows the tracked filament indicated by the blue boxes, which have
the same tracking IDs as that in previous frame in Figure 7(b).

However, in addition to the fragmentation of filaments, it is possible that
the filament maybe disappear after eruption as well as the formation of a new
filament. Therefore, in the update step, we adopt the IoU ratio to compare the
tracking result of CSRT tracker with the results of the subsequent frame after
initialization step. If the filament (merged boundary boxes) of the subsequent
frame after initialization step has the largest IoU with a certain tracked
filament by the CSRT tracker, it will be set the same tracking ID; if the IoU
is empty, i.e., there is no corresponding filament, it will be set a new ID.
In this way, the update step is finished and the result will be input to CSRT
tracker again for tracking the next subsequent frame until all frames are
tracked. The red boxes in Figure 7(d) are the results after the initialization
step, which has no IDs yet. After the update step, the red boxes in Figure
7(d) and blue boxes Figure 7(c) are compared by their IoU ratios and finally
output the tracking results indicated by the yellow boxes in Figure 7(e) and
(f). The update step can effectively handle with the situation of the eruption
of a filament. We extracted a filament with ID No.12 within field-of-view of
the blue box in Figure 8(a) as an example. We can see its dynamic evolution
and eruption from Figure 8(b) to (p) at different time, and finally
disappeared after eruption in Figure 8(q).

###  5.2 Performance of tracking

Figure 9:  The performance of our tracking method. Each blue dot represents a single filament tracking accuracy. The red triangles represent the average accuracy of each group. Table 3: Performance of our tracking method. Group number | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
Filament numbers | 15 | 14 | 19 | 12 | 18 | 17 | 20 | 13 | 22 | 21
Average accuracy | 0.81 | 0.83 | 0.81 | 0.98 | 0.80 | 0.74 | 0.82 | 0.81 | 0.78 | 0.80

We selected 10 groups of CHASE observation for testing the tracking
performance. Each group is the first day of a month and has about 15 frames of
H$\alpha$ line center images with time interval about one hour. If the
filament are tracked by our tracking method in more than 3 frames, it will be
counted for the tracking accuracy testing. If the tracking ID changed but the
ground-truth is the same filament, it would be regarded as false tracking. The
testing results are plotted in Figure 9 and summarized in Table 3. The average
tracking accuracy is $81.7\%$, which confirm the good performance of our
tracking method.

##  6 discussion and conclusion

Based on the characteristics of CHASE observations, we have developed an
efficient automated method for detecting, tracking and analyzing
filaments111The code of our detection, tracking, and analysis methods are
available on GitHub (https://github.com/ZZsolar/filament-detection-and-
tracking.git and the data are available in Solar Science Data Center of
Nanjing University (https://ssdc.nju.edu.cn/NdchaseSatellite). The code and
data have also been deposited to Zenodo under a Creative Commons Attribution
license: https://doi.org/10.5281/zenodo.10598419 (catalog
doi:10.5281/zenodo.10598419).. Instead of manually annotating filaments, we
use the K-Means method as our spectral classification tool combined with
morphological open operation to obtain the labeled filament, which guarantees
the consistence and accuracy of labelling. Then these labeled data are adopt
to train the U-Net model, and the good performance demonstrates that our
method is a viable strategy for solar filament detection. It is note that the
K-means method can also be adopted for filament detection. We adopt U-net
model for filament detection instead of K-means method due to two main
reasons. First, the K-means method is very sensitive to the uneven radiation
intensity distribution across the solar disk, which is adversely affecting
spectral classification for labeling filaments. It means that we have to
remove the limb darkening firstly. Further more, for each file, we should
select the class of the filament manually, i.e., the K-means algorithm need
additional manual assistant. So if we use K-means method for detecting, our
whole detection module would be semi-automated. Second, compared to using
H$\alpha$ spectral data as the input of K-means method during the
preprocessing steps, we only adopt the original H$\alpha$ line center image as
the input of U-Net model. In our test, the data processing and spectral
classification of a H$\alpha$ spectral file by K-means method takes several
minutes, while the detection by the U-Net model takes less than one second. It
greatly minimizes the time consumption of filament detection.

After the filament recognition, besides getting the ordinary morphological
features such as filament area and location, we apply the graph theory to
extract the filament spine, and use the cloud model to inverse the LOS
velocity distribution of the filament region. In a follow-up work, we will
implement our filament feature extraction and analysis methods to do a
statistical study on filament features, not only on their morphological
features, but also on the dynamic evolution, which are valuable for a better
understanding of the physical mechanisms of filaments. We also integrate the
U-Net model and the CSRT algorithm to track the evolution of filaments, the
test results show that our tracking method can track filaments efficiently
regardless of whether they are experiencing motions, deformations, breaking,
and eruptions.

Although our method is performed well, we also found some limitations in our
work during the experiment. First, sometimes the K-Means spectral
classification cannot figure out relative small filaments which are surrounded
by plages, as indicated in Figure 2(c) and (d). These filaments can not be
separated using intensity thresholds in H$\alpha$ line center images. The
reason may be that the K-Means algorithm only considers the distances between
different classes, while it does not account for the weights of factors such
as line depth and line width. Some other unsupervised pre-classification
models would be considered as the pre-labeling method to take account of more
information of H$\alpha$ spectra. Second, our detection method based on the
U-Net architecture, which has an excellent performance in semantic
segmentation, cannot handle the problem of filament fragments very well (if
the fragments belong to one filament without connection with each other). This
problem can affect the tracking accuracy since the filaments detected by the
U-Net model are used as inputs of the CSRT tracker. For example, if the first
frame misidentifies the filament with several fragments, the following
tracking procedure can only treat the fragments as separate filaments. The
possible ways to solve this problem require to add more information in the
labeled data and one more neural networks to merge the fragments belonging to
one filament, or adopt the instance-segmentation model like that of Guo et al.
(2022) to solve the fragment problems.

In summary, we have developed an efficient automated method for filament
detection and tracking. We utilized the U-Net model to identify filaments and
implemented the Channel and Spatial Reliability Tracking (CSRT) algorithm for
automated filament tracking. In addition, we used the cloud model to invert
the LOS velocity of filament and employed the graph theory algorithm to
extract the filament spine, which can promote understanding the dynamics of
filaments. The test favorable performance confirms the validity of our method,
which will be implemented in the following statistical analyses of filament
features and dynamics based on CHASE observations.

CHASE mission was supported by China National Space Administration. This work
was supported by NSFC under grants 12173019, 12333009, 12127901, and the CNSA
project D050101, as well as the AI & AI for Science Project of Nanjing
University.

## References

  * Asensio Ramos et al. (2023) Asensio Ramos, A., Cheung, M. C. M., Chifu, I., & Gafeira, R. 2023, Living Reviews in Solar Physics, 20, 4, doi: 10.1007/s41116-023-00038-x
  * Beckers (1964) Beckers, J. M. 1964, PhD thesis, National Solar Observatory, Sunspot New Mexico
  * Bernasconi et al. (2005) Bernasconi, P. N., Rust, D. M., & Hakim, D. 2005, Solar Phys, 228, 97, doi: 10.1007/s11207-005-2766-y
  * Bonnin et al. (2013) Bonnin, X., Aboudarham, J., Fuller, N., Csillaghy, A., & Bentley, R. 2013, Sol. Phys., 283, 49, doi: 10.1007/s11207-012-9985-9
  * Bradski (2000) Bradski, G. 2000, Dr. Dobb’s Journal of Software Tools
  * Chatzistergos, Theodosios et al. (2023) Chatzistergos, Theodosios, Ermolli, Ilaria, Banerjee, Dipankar, et al. 2023, A&A, 680, A15, doi: 10.1051/0004-6361/202347536
  * Chen (2011) Chen, P. F. 2011, Living Reviews in Solar Physics, 8, 1, doi: 10.12942/lrsp-2011-1
  * Chen et al. (2020) Chen, P.-F., Xu, A.-A., & Ding, M.-D. 2020, Research in Astronomy and Astrophysics, 20, 166, doi: 10.1088/1674-4527/20/10/166
  * Farhodov et al. (2019) Farhodov, X., Kwon, O.-H., Kang, K. W., Lee, S.-H., & Kwon, K.-R. 2019, in 2019 International Conference on Information Science and Communications Technologies (ICISCT), 1–3, doi: 10.1109/ICISCT47635.2019.9012043
  * Fuller et al. (2005) Fuller, N., Aboudarham, J., & Bentley, R. D. 2005, Solar Phys, 227, 61, doi: 10.1007/s11207-005-8364-1
  * Gao et al. (2002) Gao, J., Wang, H., & Zhou, M. 2002, Solar Phys, 205, 93
  * Gopalswamy et al. (2003) Gopalswamy, N., Shimojo, M., Lu, W., et al. 2003, ApJ, 586, 562, doi: 10.1086/367614
  * Guo et al. (2022) Guo, X., Yang, Y., Feng, S., et al. 2022, Sol. Phys., 297, 104, doi: 10.1007/s11207-022-02019-z
  * Hao et al. (2015) Hao, Q., Fang, C., Cao, W., & Chen, P. F. 2015, ApJS, 221, 33, doi: 10.1088/0067-0049/221/2/33
  * Hao et al. (2013) Hao, Q., Fang, C., & Chen, P. F. 2013, Sol. Phys., 286, 385, doi: 10.1007/s11207-013-0285-9
  * Hao et al. (2016) Hao, Q., Guo, Y., Fang, C., Chen, P.-F., & Cao, W.-D. 2016, Research in Astronomy and Astrophysics, 16, 1, doi: 10.1088/1674-4527/16/1/001
  * Hou et al. (2023) Hou, Y., Li, C., Li, T., et al. 2023, ApJ, 959, 69, doi: 10.3847/1538-4357/ad08bd
  * Kingma & Ba (2014) Kingma, D. P., & Ba, J. 2014, arXiv e-prints, arXiv:1412.6980, doi: 10.48550/arXiv.1412.6980
  * Krizhevsky et al. (2012) Krizhevsky, A., Sutskever, I., & Hinton, G. E. 2012, Advances in neural information processing systems, 25, 1097
  * Labrosse et al. (2010) Labrosse, N., Heinzel, P., Vial, J. C., et al. 2010, Space Sci. Rev., 151, 243, doi: 10.1007/s11214-010-9630-6
  * Li et al. (2019) Li, C., Fang, C., Li, Z., et al. 2019, Research in Astronomy and Astrophysics, 19, 165, doi: 10.1088/1674-4527/19/11/165
  * Li et al. (2022) —. 2022, Science China Physics, Mechanics, and Astronomy, 65, 289602, doi: 10.1007/s11433-022-1893-3
  * Lin et al. (2017) Lin, T.-Y., Goyal, P., Girshick, R., He, K., & Dollár, P. 2017, in 2017 IEEE International Conference on Computer Vision (ICCV), 2999–3007, doi: 10.1109/ICCV.2017.324
  * Liu et al. (2021) Liu, D., Song, W., Lin, G., & Wang, H. 2021, Sol. Phys., 296, 176, doi: 10.1007/s11207-021-01920-3
  * Liu et al. (2022) Liu, Q., Tao, H., Chen, C., et al. 2022, Science China Physics, Mechanics, and Astronomy, 65, 289605, doi: 10.1007/s11433-022-1917-1
  * Liu et al. (2012) Liu, R., Kliem, B., Török, T., et al. 2012, ApJ, 756, 59, doi: 10.1088/0004-637X/756/1/59
  * Lukežic et al. (2017) Lukežic, A., Vojír, T., Zajc, L. C., Matas, J., & Kristan, M. 2017, in 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 4847–4856, doi: 10.1109/CVPR.2017.515
  * MacQueen (1967) MacQueen, J. 1967, in Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume I-Theory of Statistics (University of California Press), 281–297.  https://digitalassets.lib.berkeley.edu/math/ucb/text/math_s5_v1_frontmatter.pdf
  * Martin (1998) Martin, S. F. 1998, Sol. Phys., 182, 107, doi: 10.1023/A:1005026814076
  * Mikołajczyk & Grochowski (2018) Mikołajczyk, A., & Grochowski, M. 2018, in 2018 International Interdisciplinary PhD Workshop (IIPhDW), 117–122, doi: 10.1109/IIPHDW.2018.8388338
  * Panos et al. (2018) Panos, B., Kleint, L., Huwyler, C., et al. 2018, ApJ, 861, 62, doi: 10.3847/1538-4357/aac779
  * Pierce & Slaughter (1977) Pierce, A. K., & Slaughter, C. D. 1977, Sol. Phys., 51, 25, doi: 10.1007/BF00240442
  * Qiu et al. (2022) Qiu, Y., Rao, S., Li, C., et al. 2022, Science China Physics, Mechanics, and Astronomy, 65, 289603, doi: 10.1007/s11433-022-1900-5
  * Qu et al. (2005) Qu, M., Shih, F. Y., Jing, J., & Wang, H. 2005, Solar Phys, 228, 119, doi: 10.1007/s11207-005-5780-1
  * Rezatofighi et al. (2019) Rezatofighi, H., Tsoi, N., Gwak, J., et al. 2019, Generalized Intersection over Union: A Metric and A Loss for Bounding Box Regression.  https://arxiv.org/abs/1902.09630
  * Ronneberger et al. (2015) Ronneberger, O., Fischer, P., & Brox, T. 2015, arXiv e-prints, arXiv:1505.04597, doi: 10.48550/arXiv.1505.04597
  * Shen et al. (2012) Shen, Y., Liu, Y., & Su, J. 2012, ApJ, 750, 12, doi: 10.1088/0004-637X/750/1/12
  * Shih & Kowalski (2003) Shih, F. Y., & Kowalski, A. J. 2003, Sol. Phys., 218, 99, doi: 10.1023/B:SOLA.0000013052.34180.58
  * Smith & Geach (2023) Smith, M. J., & Geach, J. E. 2023, Royal Society Open Science, 10, 221454, doi: 10.1098/rsos.221454
  * Sun et al. (2023) Sun, Z., Li, T., Tian, H., et al. 2023, ApJ, 953, 148, doi: 10.3847/1538-4357/ace5b1
  * Tian et al. (2020) Tian, Z., Shen, C., & Chen, H. 2020, in Computer Vision – ECCV 2020, ed. A. Vedaldi, H. Bischof, T. Brox, & J.-M. Frahm (Cham: Springer International Publishing), 282–298
  * Tian et al. (2023) Tian, Z., Zhang, B., Chen, H., & Shen, C. 2023, IEEE Transactions on Pattern Analysis and Machine Intelligence, 45, 669, doi: 10.1109/TPAMI.2022.3145407
  * Vial & Engvold (2015) Vial, J. C., & Engvold, O., eds. 2015, Astrophysics and Space Science Library, Vol. 415, Solar Prominences, doi: 10.1007/978-3-319-10416-4
  * Viticchié & Sánchez Almeida (2011) Viticchié, B., & Sánchez Almeida, J. 2011, A&A, 530, A14, doi: 10.1051/0004-6361/201016096
  * Wang et al. (2010) Wang, Y., Cao, H., Chen, J., et al. 2010, ApJ, 717, 973, doi: 10.1088/0004-637X/717/2/973
  * Yu et al. (2020) Yu, K., Li, Y., Ding, M. D., et al. 2020, ApJ, 896, 154, doi: 10.3847/1538-4357/ab9014
  * Yuan et al. (2011) Yuan, Y., Shih, F. Y., Jing, J., Wang, H., & Chae, J. 2011, Solar Phys, 272, 101, doi: 10.1007/s11207-011-9798-2
  * Zhang et al. (2023) Zhang, T., Hao, Q., & Chen, P. F. 2023, ApJS
  * Zhu et al. (2019) Zhu, G., Lin, G., Wang, D., Liu, S., & Yang, X. 2019, Sol. Phys., 294, 117, doi: 10.1007/s11207-019-1517-4