\documentclass[10pt]{article}
\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}
\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
linkcolor = blue,
urlcolor = blue,
citecolor = blue,
anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
\errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother
\usepackage{natbib}
\renewenvironment{abstract}
{{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
{\bigskip}
\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}
\usepackage{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}
\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}
\usepackage{float}
\begin{document}
\title{Sampling Error of Continuous Periodic Data and its Application for
Geodesy}
\author[1]{Lorant Foldvary}%
\affil[1]{Obuda University}%
\vspace{-1em}
\date{\today}
\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup
\selectlanguage{english}
\begin{abstract}
Data acquisition for geoinformatics cannot be done continuously, but by
discrete sampling of the object or phenomenon. The sampling involves
errors on the knowledge of the continuous signal due to the loss of
information in the sampling procedure. In the present study, an
analytical formulation of the sampling error is provided, which embodies
the amplitude, phase, bias and periodicity of the sampling error. The
analysis is then subsequently applied for case studies: for the GRACE
and GRACE-FO monthly solutions, and for different realizations of the
Hungarian Gravimetric Network.%
\end{abstract}%
\sloppy
\par\null
\textbf{Introduction}
Data acquisition for geoinformatics contains collection of spatially
and/or temporally changing features, which are continuous by nature. No
data acquisition can be done continuously, but by discrete sampling of
the phenomenon. Discrete sampling is often then considered to be
uninterrupted sequences of continuous data, which results in an
underestimation of the actual signal.
Shannon's proof (Shannon 1949, 1998) on the sampling theorem illustrates
elegantly the feasibility of recovering: while sampling of a
band-limited signal means a multiplication of the continuous signal with
Dirac impulses in the time domain, in the frequency domain (due to the
corresponding convolution) it yields repetition of the spectrum of the
original function. Therefore, theoretically the original function can be
recovered perfectly by filtering the sampled signals spectrum to the
original bands, which means a multiplication in the frequency domain
with a boxcar window of the proper size (see e.g. Chapter 3 of Marks
1991). This approach is, however, too idealized, as it assumes a band
limited signal, while real signals are never exactly bandlimited, but
there always an aliasing is observed. Also, in practice no ideal
anti-aliasing low-pass filters can be constructed (Unser 2000). (For
details on related issues of the sampling theory, the reader is referred
to Marks (1991, 1993) and Unser (2000)).
In fact, in practice Shannon's formulation is rarely used, as sinc
function (the time domain equivalent of the boxcar filtering) decays
slowly. For the practice, in case of an `appropriately' sampled signal,
intermediate values are assumed to be determined with `appropriate'
accuracy by linear interpolation. (Accordingly, the present study also
assumes that the sampled signal is approximated by linear
interpolation). Still, the relevance of Shannon's theorem is essential,
as (theoretically) any signal can be interpreted as infinite number of
periodic signals as Fourier transformation provides such a
decomposition.
\(f(t)=\sum_{c=0}^{\infty}{A_{c}sin(2\pi f_{c}t+\phi_{c})}\) (1)
Certainly, in practice no full equivalence of the original and the
Fourier transform signals can be achieved as the transformation can make
use of only finite numbers of frequencies, and also due to numerical
limitations. Following, however, the idealization of Shannon, in this
study, the effect of sampling is discussed for a discretization
procedure with no error, and assuming no observation errors as well.
Also, implicitly we assume the feasibility of (1), thus the discussion
focus on periodic signals only. Note also, that in the practice of
geoinformatics and Remote Sensing the discretization is often obtained
by averaging over a finite segment of the data, e.g. Digital Terrain
Models (DTM) are determined based on several observations referring to a
certain pixel, resulting in aliasing the point-wise data by the block
averaging (F\selectlanguage{ngerman}öldváry 2015). This study assumes perfectly point-wisely
sampled data.
Nowadays, sampling algorithms are generalizing the Shannonian approach:
instead of constraining the signal into a limited bandwidth before
recovering it, the proper band-limited filter can be achieved by methods
applying orthogonal, frequency dependent base functions, such as
wavelets (Strang-Nguyen 1996; Mallat 1998), finite elements (Strang
1971; Selesnick 1999), frames (Duffin-Schaeffer 1952; Benedetto 1992),
among others.
In summary, the aim of the study is to provide a mathematical tool for
sampling error estimation. The sampling, by its discrete characteristics
involves errors on the knowledge of the real continuous phenomenon.
Without understanding the limitations of the discretization, the
observed phenomenon may be interpreted falsely. In the present study,
analytical formulation of the sampling error is to be provided, which
embodies the characteristics of the sampling error by determining its
amplitude, phase, bias and periodicity. Such information can be used for
planning optimal resolution of sampling a process but cannot be used
(and it is not the scope of this study) for reconstructing the original
signal.
\textbf{Formulation}
Figure 1 displays an example of a periodic signal, which is sampled with
a certain sampling period. Let \(T\) refer to the period of
the signal, and \(T\) to the sampling period. If the
sampling period is sufficiently fine, then the signal can efficiently be
approximated by assuming linearity between two consecutive points. In
practice, for sake of simplicity, such a linearity is assumed;
accordingly, linear interpolation is used for approximating inner values
of the signal.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}
\end{center}
\end{figure}
\textbf{Figure 1.} Sampled signal is approximated by linear
interpolation. (The periodic signal is represented by the blue curve,
the sampled values are red circles, and the linearly interpolated signal
is the red dashed line).
Assume that the periodic signal with a period of \(T\)
(equivalently can be described by its frequency, \(\omega=1/T\) ) is
sampled regularly by a constant sampling rate, \(T=1/f_{s}\), where
\(f_{s}\) is the sampling frequency. The dimension of the
signal may range at different scales, even the content of the signal can
be different, e.g. time in seconds, hours, years, ages or length in mm,
m, km, AU. Therefore, instead of discussing actual reliable scales and
magnitudes, in the present analysis the periodic function is
characterized by its amplitude, \(A\) and frequency,
\(\omega\), while the sampling is defined by the sampling rate,
\(T\), or equivalently by the ratio of the frequency of the
signal, \(\omega\) and the observation, \(f_{s}\)
\(N=\frac{\ \omega}{f_{s}}=\frac{T}{T}\) (2).
The amplitude, \(A\) is set arbitrarily to a unit, and
error estimation due to the sampling is performed by considering two
parameters: the frequency, \(\omega\) of the signal and the
ratio \(N\). The error estimates are provided in
percentage of the amplitude, \(A\).
Beyond\(A\), \(\omega\) and \(N\),
also the phase of the signal, \(\phi\) is used for defining
the periodic signal, however the calculus later is performed
independently of this variable. All in all, the periodic function is
defined as
\(f(x)=Asin(2\pi\omega\ x+\phi)\) (3),
where \emph{x} is the independent variable. The sampling affects the
knowledge of the observed signal between the sampling epochs. Therefore,
sampling error is modelled here as the difference of real and the
(linearly) interpolated values. According to Figure 2, when the function
value at an arbitrary epoch, \emph{x\textsubscript{k}} falling into the
interval of {[}\emph{x\textsubscript{i}} ,\emph{x\textsubscript{i+1}}
{]}, the error due to the sampling becomes the difference of the real,
\emph{f(x\textsubscript{k})} and the interpolated,
\emph{f\textsubscript{int}(x\textsubscript{k})} function values,
\(\varepsilon_{k}\).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}
\end{center}
\end{figure}
\textbf{Figure 2.} Errors due to sampling a periodic signal and
approximating by linear interpolation
The interpolated function
value,\emph{f\textsubscript{int}(x\textsubscript{k})} can be derived by
considering it as a point of division, i.e. it divides in a ratio
of\emph{(x\textsubscript{k}-x\textsubscript{i})}
:\emph{(x\textsubscript{i+1}-x\textsubscript{k})} the join of points
(\emph{x\textsubscript{i}} , \emph{f(x\textsubscript{i})} ) and
(\emph{x\textsubscript{i+1}} , \emph{f(x\textsubscript{i+1})} ):
\(f_{\text{int}}\left(x_{k}\right)=\frac{\left(x_{i+1}-x_{k}\right)f\left(x_{i}\right)+\left(x_{k}-x_{i}\right)f\left(x_{i+1}\right)}{x_{i+1}-x_{i}}\)(4).
Accordingly, the sampling error at this epoch becomes
\(\varepsilon_{k}=f\left(x_{k}\right)-f_{\text{int}}\left(x_{k}\right)=f(x_{k})-\frac{\left(x_{i+1}-x_{k}\right)f\left(x_{i}\right)+\left(x_{k}-x_{i}\right)f\left(x_{i+1}\right)}{x_{i+1}-x_{i}}\)(5).
Note that the difference of two consecutive abscissae (in the
denominator of the second term of the righthand side) is the sampling
interval, \(T=x_{i+1}-x_{i}\).
Generally, the sampling error in the engineering practice is estimated
by L2-norm; for some applications also L1-norm and L-norm is used. The
present study focuses on the L2-norm, which is derived after the simpler
form of L1-norm. The L1-norm of the sampling error over the
{[}\emph{x\textsubscript{i}} , \emph{x\textsubscript{i+1}} {]} interval
can be defined as it is the mean of the \(\varepsilon\) differences
at all epochs, i.e.
\(\sigma_{L1}\left([x_{i},x_{i+1}]\right)=\frac{L1\left([x_{i},x_{i+1}]\right)}{n}=\frac{\sum_{x_{i}}^{x_{i+1}}\left|\varepsilon_{k}\right|}{n}=\frac{\int_{x_{i}}^{x_{i+1}}{\left|\varepsilon(x)\right|\text{dx}}}{T}\)(6).
The integral provides the area between the real and interpolated curves,
which is then subsequently divided by the sampling interval. By rigorous
calculus, a closed-form for the L1-norm (inserting (3) and (5) to (6))
can be derived by solving the definite integral of
\(\sigma_{L1}\left([x_{i},x_{i+1}]\right)=\frac{\int_{x_{i}}^{x_{i+1}}{\left|Asin(2\pi\omega\ x+\phi)-\frac{\left(x_{i+1}-x\right)Asin(2\pi\omega\ x_{i}+\phi)+\left(x-x_{i}\right)Asin(2\pi\omega\ x_{i+1}+\phi)}{T}\right|\text{dx}}}{T}\)(7),
resulting in
\(\sigma_{L1}\left([x_{i},x_{i+1}]\right)=\left|C\bullet\cos\left(2\pi\omega x_{i}+\phi\right)+S\bullet sin(2\pi\omega x_{i}+\phi)\right|\)(8),
where
\(C=\frac{A}{2\pi\omega T}-\frac{A}{2}sin(2\pi\omega T)-\frac{A}{2\pi\omega T}\cos\left(2\pi\omega T\right)\)(9)
and
\(S=-\frac{A}{2}+\frac{A}{2\pi\omega T}\sin\left(2\pi\omega T\right)-\frac{A}{2}cos(2\pi\omega T)\)(10).
For details of deriving (8), see Appendix A. By introducing the
following notation
\(n=2\pi N=2\pi\omega T\) (11)
it simplifies to
\(C=\frac{A}{n}-\frac{A}{2}sin(n)-\frac{A}{n}\cos\left(n\right)\)(12)
and
\(S=-\frac{A}{2}+\frac{A}{n}\sin\left(n\right)-\frac{A}{2}cos(n)\)(13).
Using the\(C\bullet\cos\left(\alpha\right)+S\bullet\sin\left(\alpha\right)=R\bullet\sin\left(\alpha+\varphi\right)\)conversion, the parameters
\(C\) and \(S\) can be replaced to
\(R\) and\(\phi\) as
\(\sigma_{L1}\left([x_{i},x_{i+1}]\right)=R\bullet\left|sin(2\pi\omega x_{i}+\phi+\varphi)\right|\)(14),
where
\(R=\sqrt{2}\frac{A}{n}\sqrt{\left(1+\frac{n^{2}}{4}\right)-nsin(n)-\left(1-\frac{n^{2}}{4}\right)\cos\left(n\right)}\)(15)
and
\(\varphi=arctg\frac{1-\frac{n}{2}\sin\left(n\right)-cos(n)}{-\frac{n}{2}+\sin\left(n\right)-\frac{n}{2}cos(n)}\)(16).
Note that (15) and (16) depends purely on the ratio of the sampling
interval and the period of the signal, \(N\), c.f. (11).
In (15) it is also multiplied by \(A\), which provides the
scale (and the unit) of\(R\), i.e. it works as the scale
factor of the ordinate of the figure, while only the ratio
\(N\) is relevant along the abscissa.
Similarly, the sampling error based on the L2-norm can also be derived
applying the definition of RMS,
\(\sigma_{L2}\left([x_{i},x_{i+1}]\right)=\frac{L2\left([x_{i},x_{i+1}]\right)}{\sqrt{n}}=\sqrt{\frac{\sum_{x_{i}}^{x_{i+1}}{\varepsilon_{k}}^{2}}{n}}=\sqrt{\frac{\int_{x_{i}}^{x_{i+1}}{{\varepsilon(x)}^{2}\text{dx}}}{T}}\)(17).
Making use of (3) and (5), the definite integral in (17) becomes
\begin{equation}
{L2\left(\left[x_{i},x_{i+1}\right]\right)}^{2}=\nonumber \\
\end{equation}
\(\int_{x_{i}}^{x_{i+1}}{\left(Asin(2\pi\omega\ x+\phi)-\frac{\left(x_{i+1}-x\right)Asin(2\pi\omega\ x_{i}+\phi)+\left(x-x_{i}\right)Asin(2\pi\omega\ x_{i+1}+\phi)}{T}\right)^{2}\text{dx}}\)(18),
following similar steps of derivation to Appendix A, it can be expressed
in a compact form of
\({L2\left(\left[x_{i},x_{i+1}\right]\right)}^{2}=C\bullet\cos\left(4\pi\omega x_{i}+2\phi\right)+S\bullet\sin\left(4\pi\omega x_{i}+2\phi\right)+B\)(19),
where
\begin{equation}
C=-\frac{A^{2}T}{24\pi^{2}\omega^{2}{T}^{2}}\left\{4\pi^{2}\omega^{2}{T}^{2}-6+\left(4\pi^{2}\omega^{2}{T}^{2}-6\right)\cos\left(4\pi\omega T\right)+\left(4\pi^{2}\omega^{2}{T}^{2}+12\right)\cos\left(2\pi\omega T\right)-9\pi\omega Tsin(4\pi\omega T)\right\}\nonumber \\
\end{equation}
(20),
\begin{equation}
S=-\frac{A^{2}T}{24\pi^{2}\omega^{2}{T}^{2}}\left\{9\pi\omega T-\left(4\pi^{2}\omega^{2}{T}^{2}-6\right)\sin\left(4\pi\omega T\right)-\left(4\pi^{2}\omega^{2}{T}^{2}+12\right)\sin\left(2\pi\omega T\right)-9\pi\omega Tcos(4\pi\omega T)\right\}\nonumber \\
\end{equation}
(21),
and
\(B=-\frac{A^{2}T}{24\pi^{2}\omega^{2}{T}^{2}}\left\{12-20\pi^{2}\omega^{2}{T}^{2}-\left(4\pi^{2}\omega^{2}{T}^{2}+12\right)\cos\left(2\pi\omega T\right)\right\}\)(22).
By introducing the \(n\) defined by (11), it simplifies to
\begin{equation}
C=-\frac{A^{2}T}{6n^{2}}\left\{n^{2}-6+\left(n^{2}-6\right)\cos\left(2n\right)+\left(n^{2}+12\right)\cos\left(n\right)-\frac{9}{2}nsin(2n)\right\}\nonumber \\
\end{equation}
(23),
\begin{equation}
S=-\frac{A^{2}T}{6n^{2}}\left\{\frac{9}{2}n-\left(n^{2}-6\right)\sin\left(2n\right)-\left(n^{2}+12\right)\sin\left(n\right)-\frac{9}{2}ncos(2n)\right\}\nonumber \\
\end{equation}
(24),
and
\(B=-\frac{A^{2}T}{6n^{2}}\left\{12-{5n}^{2}-\left(n^{2}+12\right)\cos\left(n\right)\right\}\)(25).
Using the\(C\bullet\cos\left(\alpha\right)+S\bullet\sin\left(\alpha\right)=R\bullet\sin\left(\alpha+\varphi\right)\)conversion, the parameters
\(C\) and \(S\) can be replaced to
\(R\) and\(\varphi\) as
\({L2\left([x_{i},x_{i+1}]\right)}^{2}=R\bullet\sin\left(4\pi\omega x_{i}+2\phi+\varphi\right)+B\)(26),
where
\(R=\frac{A^{2}T}{6n^{2}}\left(\left(n^{2}+12-9n\sin\left(n\right)+\left(2n^{2}-12\right)\cos\left(n\right)\right)^{2}\operatorname{}\left(n\right)+\frac{1}{4}\left(-9n+9n\cos\left(2n\right)+\left(2n^{2}-12\right)\sin\left(2n\right)+\left(2n^{2}+24\right)\sin\left(n\right)\right)^{2}\right)^{\frac{1}{2}}\)
(27),
\(\varphi=arctg\frac{n^{2}-6+\left(n^{2}-6\right)\cos\left(2n\right)+\left(n^{2}+12\right)\cos\left(n\right)-\frac{9}{2}nsin(2n)}{\frac{9}{2}n-\left(n^{2}-6\right)\sin\left(2n\right)-\left(n^{2}+12\right)\sin\left(n\right)-\frac{9}{2}ncos(2n)}\)(28)
and \(B\) is as in (25). As (23) to (28) refers
to\({L2\left(\left[x_{i},x_{i+1}\right]\right)}^{2}\), the sampling error can be determined by (17) as
\(\sigma_{L2}\left([x_{i},x_{i+1}]\right)=\frac{L2\left([x_{i},x_{i+1}]\right)}{\sqrt{T}}=\sqrt{\frac{{L2\left(\left[x_{i},x_{i+1}\right]\right)}^{2}}{T}}=\sqrt{\frac{R\bullet\sin\left(4\pi\omega x_{i}+2\phi+\varphi\right)+B}{T}}\)(29),
so no closed-form solution can be derived. Similarly to (15) and (16),
(27) and (28) depends purely on the ratio of the sampling interval and
the period of the signal, \(N\), and the amplitude,
\(R\) provides the scale and the unit in (27) by the
multiplicator of \(A^{2}T\).
\textbf{Generalization of the formulation}
A major shortcoming of the derived formulation is that it holds for
purely periodic signal, which is indeed unrealistic. Note, however, that
theoretically it can be generalized to any time series, as a Fourier
series representation of a function consists of infinite number of
periodic components, c.f. equation (1), where summation is done
by\emph{c} , and the reciprocal of any
wavelengths,\emph{T\textsubscript{c}} , is the
frequency,\(f_{c}=\frac{1}{T_{c}}\). Due to the orthonormality of the Fourier
base functions, the effect if the sampling error can be estimated to
each component independently along the Fourier spectrum. Considering the
frequency of a Fourier component of the signal, \(f_{c}\) and
that of the sampling, \(f_{s}\), the ratio \emph{N} for each
component can be determined similarly to equation (2):
\(N_{c}=\frac{\ f_{c}}{f_{s}}=\frac{T}{T_{c}}\) (30).
It can then be used for each Fourier components independently using the
equations of section 2.
Note that the bias does not need special attention in this case. As for
the \emph{c=0} term the Fourier component \(T_{0}=\infty\)
and\(f_{0}=0\), for any non-zero sampling frequency,
\(f_{s}\neq 0\), equation (30) becomes zero. For such a case,
equations (11) to (16) and (23) to (29) has singularity, therefore the
signal should be unbiased before estimating the sampling error.
\textbf{Discussion}
When sampling error of a periodic signal is to be estimated, it should
definitely be kept in mind that it highly depends on which part of the
signal it takes place, c.f. in Figure 1 the errors are apparently larger
at the extremes (peaks) than around the inflection points.
Equations (8) to (10) and (19) to (22) provide analytic formulations for
the periodicity of the sampling error. These equations indicate that it
purely depends on the ratio of the signal period and the sampling
interval, \(N\). The periodicity of the sampling error is
displayed in Figure 3 for an arbitrary example of N=0.05 ratio, i.e. 20
samplings during a period. It shows that the periodicity of the error is
tied to the periodicity of the signal, but with a phase lag (it is 9\selectlanguage{ngerman}° in
this example). It also shows that for such a case the maximal error is
0.82\% and 0.90\% of the amplitude of the period of the signal,
respectively for the L1 and L2 norms.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}
\end{center}
\end{figure}
\textbf{Figure 3.} Sampling error due of a periodic signal with N=0.05
ratio.
The L1-norm and the L2-norm based sampling errors are shown in Figure 4
for a range of \(N\) ratios in the interval of {[}0.01,
0.5{]}. The limits were chosen considering practical aspect. The tested
finest resolution was 0.01 (equivalent to 10 samplings per each period),
which has resulted in negligible, 0.03\% and 0.04\% errors using the
L1-norm and the L2-norm, respectively. The tested largest ratio was
chosen in accordance with the Nyquist criterion, which says that
meaningful frequency components of the properly sampled periodic signal
exist below the half of the period of the signal, i.e. the Nyquist
frequency (see e.g. Leis, 2011). For this extreme case, the sampling
error was found to be 63.7\% and 70.7\% for the L1-norm and the L2-norm,
respectively.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}
\end{center}
\end{figure}
\textbf{Figure 4.} Sampling error for ratios ranges from 0.01 to 0.5.
\textbf{Example in the time domain: GRACE monthly solutions}
Temporal variations of the gravity field globally are monitored by the
GRACE (between 2002-2017) and the GRACE-FO (from 2018) missions
(Bettadpur, 2007). Mass variation in the Amazonas basin is dominated by
the annual period of the hydrological cycle (c.f. gravity anomaly time
series from GRACE and GRACE-FO monthly solutions in Figure 5). The
monthly sampling is equivalent to N=1/12=0.083. Using (14) and (29), the
maximal L1-norm sampling error becomes 2.27\%, while the maximal L2-norm
error is 2.49\%.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}
\end{center}
\end{figure}
\textbf{Figure 5.} Gravity anomaly in the Amazonas basin using GRACE and
GRACE-FO monthly solutions.
According to F\selectlanguage{ngerman}öldváry (2015), the annual amplitude is underestimated due
to the averaging over the sampling period by 1.15\%, thus the effect of
the sampling error exceeds this phenomenon. Nevertheless, the GRACE and
GRACE-FO models notably scatter from the annual characteristics (c.f.
Figure 5, where also a best fit annual and semiannual curve is
presented). The deviations of the detrended and unbiased GRACE/GRACE-FO
monthly gravity anomalies from an annual and semiannual periodic curve
is 35.58\% of the magnitude of the signal, in average. This deviation
consists of non-periodic signal, modelling errors (such as errors of the
atmospheric correction, leakage of signal from outside of the basin,
errors of smoothing and de-correlation filtering) and observation errors
(integrating such errors as system-noise error in the KBR range-rate
measurements, accelerometer error, errors of the ultra-stable
oscillator, and orbit errors) as well (Swenson et al., 2003).
Considering that the errors of GRACE are believed to be largely free of
annual components (Wahr et al., 2006), the estimated annual component on
Figure 5 can be considered to be purely due to the annual variations of
the total water storage, consequently, it can be considered to be a real
signal, while the error content is included in the non-periodic terms.
Accordingly, the unmodelled non-periodic components and observation
errors can be handled independently from the annual component, and the
later component can be described by an error of 2.49\% of the signal
content due to the sampling and an error of 1.15\% of the signal content
due to the averaging.
\textbf{Example in the space domain: the Hungarian Gravimetric Network}
In order to indicate the relevance of the errors contaminated by the
sampling, a case study of the comparison of two epochs of the Hungarian
Gravimetric Network (MGH) is provided. The first epoch of the network
was dated to the 1950s, labelled as MGH-50 (Facsinay-Szilárd, 1956),
while the second one is dated to 2000, referred to as MGH-2000 (Csapó,
2000). The distribution of the gravimetric stations is shown in Figure 6
for the MGH-50 and the MGH-2000 networks.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}
\end{center}
\end{figure}
\textbf{Figure 6.} The point distribution of the MGH-50 and MGH-2000
gravimetric network
The covered area (the area of the country) is approximately 93000
km\textsuperscript{2}, the number of stations is more than 900 for both
the MGH-50 and the MGH-2000 networks, the gravity interval of the
country is approximately 350 mGal. The distance between gravimetric
stations is roughly 10 km, thus it provides a nice coverage of the
country for capturing the relevant features of the gravity field.
For historical reasons, however, most stations of the MGH-50 network
have been perished, accordingly, most points of the MGH-2000 are newly
established stations. In fact, there are only 13 identical stations of
the two networks (see common points in Figure 6). Apart from the change
of the gravity by time (which is not striking due to the moderate
tectonics of the area), the two networks are meant to present two
independent description of the same phenomenon, that is the gravity
field of Hungary.
Based on the gravity values of MGH-50 and MGH-2000 networks, gravity has
been interpolated onto a common grid, and the difference of the two
grids is displayed in Figure 7. Even though the difference between them
is in the range of some mGals, there are striking outliers reaching even
the 100 mGal value. This is very much, which cannot be interpreted
neither with the temporal variations of the gravity field nor with
observation errors. This difference is obviously a consequence of the
different sampling of the gravity field.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}
\end{center}
\end{figure}
\textbf{Figure 7.} Differences of interpolated gravity fields based on
MGH-50 and on MGH-2000 gravimetric network data. The values of the
colourbar are in mGal unit.
Using the formulations above, the adequacy of the resolution of the
Hungarian gravity network can be tested. According to equation (15),
with stations located about 10 km distance from each other, features of
the gravity field can be observed with 5\%, 1\% and 0.1\% sampling error
up to a resolution of 0.18 km, 0.42 km and 1.33 km, respectively.
Indeed, local gravity anomalies can be relevant with even finer
resolutions, therefore it is obvious that the 10 km distance between the
stations are capable only for capturing middle wavelength gravity field
variations.
A more reliable and informative test can be provided if also the
magnitude of the observable is considered. Since the gravity field sums
up of gravity features at different scales (ranging from local to
global), the spectral behaviour of the gravity field should be modelled.
In our case it is approximated by the Kaula's rule of thumb (Kaula
1966), which is expressed in equation (31) for gravity anomaly,
\(g\):
\(\sigma\left(g\right)=\frac{\text{kM}}{R^{2}}\left(n-1\right)\bullet\frac{\sqrt{2n+1}}{n^{2}}10^{-5}\)(31).
The Kaula's rule of thumb estimates the signal content of gravity as a
function of the spherical harmonic degree, \(n\), for an
Earth with an average mass \(M\) and an average radius
\(R\). The `classical' rule of thumb (referring to the
potential) is converted to gravity anomaly,\(g\) by
multiplying it with the transfer function of the spherical harmonic
expansion of the gravity anomaly. (For detailed explanation on the
physical content of the Kaula's rule of thumb see Kaula (1992), Rummel
(2004) and McMahon et al. (2016)).
The spherical harmonic degree, \(n\) can approximately be
related to the scale of the gravity content, i.e. to the spatial
resolution of the features of the gravity field (described by a
wavelength \(T\)).
\(T\left(g\right)=\frac{2\pi R}{n}\) (32).
Using (31) and (32), Figure 8 shows an estimate of the gravity field
content by the spatial resolution according to Kaula's rule of thumb in
logarithmic scale. In the discussion below, two terms highly used in
geodesy are applied here. These are: omission error refers to those
errors, which are committed due to omitting certain parts of the
frequency spectrum, while commission error refers to the errors provided
by the involved frequencies. Certainly, the total error is a combination
of the omission and commission errors, and one can only attempt to
separate them based on the spatial resolution of the observed total
error.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image8/image8}
\end{center}
\end{figure}
\textbf{Figure 8.} Signal content of gravity anomaly based on Kaula's
rule of thumb vs. sampling error (by L1-norm and L2-norm as well)
according to a sampling rate, \(T\) of 10 km.
The sampling error is estimated by L1-norm and L2-norm using equation
(15). The value of the sampling rate has been set to \(T=10km\),
and the values of \(N\) were determined by using
\(T\) as the independent variable. The resulted curves are
shown in Figure 8. The maximal value of \(N\) was 0.5
according to the Nyquist criterion, which is at the distance of 20 km.
The \(N=1\) abscissa is displayed on the figure by a vertical
dotted line, which refers to the \(T=T=10km\) value. By any mean,
signal content below \(T=10km\) appears as noise in the signal,
as the corresponding frequencies are overlooked (not sampled) by the
\(T=10km\) sampling rate. The omission of the signal content of
the short wavelength gravity (i.e. at less than \(10km\)
resolution) may particularly be relevant at mountainous regions, where
the gravity may change more sharply, the gradients of the gravity are
larger. Basically, omission error of the short wavelength gravity is the
primer source of 100 mGal large, very localized outliers in Figure 7.
Beyond the omission error, the commission error due to the sampling is
also large: according to Figure 8, the sampling generated error reaches
the signal at 28.1 km (L1-norm) and 26.8 km (L2-norm). A signal is
reliable up to that point where its noise content is at least one order
of magnitude smaller than the signal content, i.e. the Signal-to-Noise
ratio is 10. In the case of Figure 8, this is reached at 73.8 km
(L1-norm) and 70.7 km (L2-norm), so any smaller scale gravity
information can be gained only quite uncertainly due to the commission
error. As in an actual case both the omission and commission errors
affect the solution, the gravity field with 10 km sampling may describe
properly the gravity field features up to 100 km resolution only.
\textbf{Summary}
Sampling of continuous signals cannot be avoided for practical
applications therefore sampling errors are contaminating the knowledge
of the continuous signal. In order to check whether the sampling is
sufficiently fine for a certain application, a simple test on the
sampling rate can be performed. In this study, analytic formulation for
L1-norm sampling error estimate of a periodic signal has been delivered
in a closed form. Also, the L2-norm error estimate has been derived
making use of a symbolic programming module of Matlab, which has not
been presented here but applied for the tests. According to the results,
the 5\%, 1\% and 0.1\% errors can be reached by N=0.177, 0.078 and
0.025, equivalent to 6, 13 and 41 samples per period, respectively.
As an example, GRACE-borne gravity anomaly time series are analysed for
the region of the Amazonas basin. It was found that the unmodelled
non-periodic components and observation errors can be handled
independently from the annual component, and that the annual component
can be described by an error of 2.49\% of the signal content due to the
sampling and an error of 1.15\% of the signal content due to the
averaging.
As a spatial example, the case of the Hungarian gravity network is
analysed. With stations located about 10 km distance from each other,
both omission and commission errors were found to be relevant. Fine
features of the gravity field are largely contaminated by the omission
of the gravity signal over less than 10 km scales, resulting in huge,
local peaks in the gravity anomaly error map in Figure 6. Commission
error was found to be relevant up to a resolution of approximately 70
km. As in an actual case both the omission and commission errors affect
the solution, the gravity field with 10 km sampling may describe
properly the gravity field features up to 100 km resolution only. For
reducing the spatial resolution, the sampling density should be
increased further.
\textbf{References}
\begin{enumerate}
\tightlist
\item
Benedetto, J. J. (1992): Irregular sampling and frames, in: Wavelets -
A Tutorial in Theory and Applications, C.K. Chui (Ed.), Boca Raton,
FL, CRC Press, pp. 445--507.
\item
Bettadpur, S. (2007): Gravity Recovery and Climate Experiment Level-2
Gravity Field Product User Handbook, GRACE 327-734, CSR Publ.
GR-03-01, Rev 2.3, pp. 19, University of Texas at Austin, USA.
\item
Csap\selectlanguage{ngerman}ó, G. (2000): The Hungarian Gravimetric Network, MGH-2000 (In
Hungarian; Magyarország új gravimetriai alaphálózata, MGH--2000),
Geodézia és Kartográfia, 52(2), p. 27-33.
\item
Duffin, R. J., Schaeffer, A. C. (1992): A class of nonharmonic Fourier
series, Trans. Amer. Math. Soc., vol. 72, pp. 314--366.
\item
Facsinay L., Szilárd J. (1956): The Hungarian Gravimetric Network (in
Hungarian; A magyar országos gravitációs hálózat), Geofizikai
Közlemények, V(2), p. 3-49.
\item
Földváry, L (2015).: Desmoothing of averaged periodical signals for
geodetic applications, Geophysical Journal International, 201 (3):
1235-1250, DOI 10.1093/gji/ggv092
\item
Kaula, W, M. (1966): Theory of Satellite Geodesy, Blaisdell, Waltham
\item
Kaula W.M. (1992): Properties of the Gravity Fields of Terrestrial
Planets. In: Colombo O.L. (ed.) From Mars to Greenland: Charting
Gravity With Space and Airborne Instruments. International Association
of Geodesy Symposia, vol 110. Springer, New York, NY, 1-10.
\item
Leis, J. W. (2011): Digital Signal Processing Using MATLAB for
Students and Researchers. John Wiley \& Sons. p. 82, ISBN
9781118033807
\item
Mallat, S. (1998): A Wavelet Tour of Signal Processing. San Diego, CA,
Academic.
\item
Marks, R.J. (1991): Introduction to Shannon Sampling and Interpolation
Theory, Springer-Verlag
\item
Marks, R.J., editor (1993): Advanced Topics in Shannon Sampling and
Interpolation Theory, Springer-Verlag
\item
McMahon, J. W., Farnocchia, D., Scheeres, D., Chesley, S. (2016):
Understanding Kaula's Rule for Small Bodies, 47th Lunar and Planetary
Science Conference 2016, paper 2129.
\item
Rummel, R. (2004): Gravity and Topography of Moon and Planets. Earth
Moon Planet 94, 103--111. https://doi.org/10.1007/s11038-005-3245-z
\item
Selesnick, I. W. (1999): Interpolating multiwavelets bases and the
sampling theorem, IEEE Trans. Signal Processing, vol. 47, no. 6, pp.
1615--1621.
\item
Shannon, C. E. (1949): Communication in the presence of noise, in:
Proc. IRE, vol. 37, pp. 10-21.
\item
Shannon, C. E. (1998): Classic paper: Communication in the presence of
noise, Proc. IEEE, vol. 86, no. 2, pp. 447--457.
\item
Strang, G. (1971): The finite element method and approximation theory,
in: Numerical Solution of Partial Differential Equations---II, B.
Hubbard (Ed.), New York, Academic, pp. 547--583.
\item
Strang, G., Nguyen, T. (1996): Wavelets and Filter Banks. Wellesley,
MA, Wellesley-Cambridge.
\item
Swenson, S., J. Wahr, and P. C. D. Milly (2003): Estimated accuracies
of regional water storage variations inferred from the gravity
recovery and climate experiment (GRACE), Water Resources Research, 39,
11--1. https://doi.org/10.1029/2002WR001808
\item
Unser, M. (2000): Sampling-50 Years after Shannon, Proc. IEEE, 88(4),
p. 569--587.
\item
Wahr, J., S. Swenson, and I. Velicogna (2006), Accuracy of GRACE mass
estimates, Geophysical Research Letteres, 33, L06,401,
doi:10.1029/2005GL025305.
\end{enumerate}
\textbf{Appendix A}
According to equation (8), a closed-form for the L1-norm can be derived
by solving equation (7), which reads
\(\sigma_{L1}\left([x_{i},x_{i+1}]\right)=\frac{\int_{x_{i}}^{x_{i+1}}{\left|Asin(2\pi\omega\ x+\phi)-\frac{\left(x_{i+1}-x\right)Asin(2\pi\omega\ x_{i}+\phi)+\left(x-x_{i}\right)Asin(2\pi\omega\ x_{i+1}+\phi)}{T}\right|\text{dx}}}{T}\)(A1).
Concentrating on the integration in the numerator, there are two
integrations should be completed, since
\(L1\left(\left|\varepsilon\right|\right)=L1\left(\varepsilon\right)\ \text{or}\ L1\left(\left|\varepsilon\right|\right)=L1\left(-\varepsilon\right)\ \)(A2).
In the\(L1\left(\left|\varepsilon\right|\right)=L1\left(\varepsilon\right)\)case the derivation starts with splitting the
integrand into separated terms,
\(\int_{x_{i}}^{x_{i+1}}{\left(Asin(2\pi\omega\ x+\phi)-\frac{\left(x_{i+1}-x\right)Asin(2\pi\omega\ x_{i}+\phi)+\left(x-x_{i}\right)Asin(2\pi\omega\ x_{i+1}+\phi)}{T}\right)\text{dx}}=\int_{x_{i}}^{x_{i+1}}{A\sin\left(2\pi\omega\ x+\phi\right)\text{dx}}-\int_{x_{i}}^{x_{i+1}}{\frac{x_{i+1}\operatorname{Asin}\left(2\pi\omega\ x_{i}+\phi\right)}{T}\text{dx}}+\int_{x_{i}}^{x_{i+1}}{\frac{\operatorname{Asin}\left(2\pi\omega\ x_{i}+\phi\right)}{T}\text{xdx}}-\int_{x_{i}}^{x_{i+1}}{\frac{\operatorname{Asin}\left(2\pi\omega\ x_{i+1}+\phi\right)}{T}\text{xdx}}+\int_{x_{i}}^{x_{i+1}}{\frac{x_{i}\operatorname{Asin}\left(2\pi\omega\ x_{i+1}+\phi\right)}{T}\text{dx}}\)(A3).
By performing the definite integration, and making use of
\(x_{i+1}=x_{i}+\text{T\ }\) (A4),
equation (A3) becomes
\(\left[-\frac{A}{2\pi\omega}\cos\left(2\pi\omega\ x+\phi\right)\right]_{x_{i}}^{x_{i+1}}-\left[\frac{x_{i+1}\operatorname{Asin}\left(2\pi\omega\ x_{i}+\phi\right)}{T}x\right]_{x_{i}}^{x_{i+1}}+\left[\frac{\operatorname{Asin}\left(2\pi\omega\ x_{i}+\phi\right)}{T}\frac{x^{2}}{2}\right]_{x_{i}}^{x_{i+1}}-\left[\frac{\operatorname{Asin}\left(2\pi\omega\ x_{i+1}+\phi\right)}{T}\frac{x^{2}}{2}\right]_{x_{i}}^{x_{i+1}}+\left[\frac{x_{i}\operatorname{Asin}\left(2\pi\omega\ x_{i+1}+\phi\right)}{T}x\right]_{x_{i}}^{x_{i+1}}=-\left(\frac{A}{2\pi\omega}\cos\left(2\pi\omega\ x_{i}+2\pi\omega\ T+\phi\right)-\frac{A}{2\pi\omega}\cos\left(2\pi\omega\ x_{i}+\phi\right)\right)-\left(x_{i}+T\right)\operatorname{Asin}\left(2\pi\omega\ x_{i}+\phi\right)+\operatorname{Asin}\left(2\pi\omega\ x_{i}+\phi\right)\left(x_{i}+\frac{T}{2}\right)-\operatorname{Asin}\left(2\pi\omega\ x_{i}+2\pi\omega\ T+\phi\right)\left(x_{i}+\frac{T}{2}\right)+x_{i}\operatorname{Asin}\left(2\pi\omega\ x_{i}+2\pi\omega\ T+\phi\right)=-\frac{A}{2\pi\omega}\cos\left(2\pi\omega\ x_{i}+2\pi\omega\ T+\phi\right)+\frac{A}{2\pi\omega}\cos\left(2\pi\omega\ x_{i}+\phi\right)-\frac{T}{2}\operatorname{Asin}\left(2\pi\omega\ x_{i}+\phi\right)-\frac{T}{2}\operatorname{Asin}\left(2\pi\omega\ x_{i}+2\pi\omega\ T+\phi\right)\)(A5).
Applying trigonometric identities
\(\cos\left(\alpha+\beta\right)=\cos\alpha\cos\beta-\sin\alpha\sin\selectlanguage{greek}\text{β\ }\selectlanguage{english}\)(A6)
and
\(\sin\left(\alpha+\beta\right)=\sin\alpha\cos\beta+\cos\alpha\sin\selectlanguage{greek}\text{β\ }\selectlanguage{english}\)(A7),
equation (A5) is reformulated as
\(-\frac{A}{2\pi\omega}\left(\cos\left(2\pi\omega\ x_{i}+\phi\right)\cos\left(2\pi\omega\ T\right)-\sin\left(2\pi\omega\ x_{i}+\phi\right)\sin\left(2\pi\omega\ T\right)\right)+\frac{A}{2\pi\omega}\cos\left(2\pi\omega\ x_{i}+\phi\right)-\frac{T}{2}\operatorname{Asin}\left(2\pi\omega\ x_{i}+\phi\right)-\frac{T}{2}A\left(\sin\left(2\pi\omega\ x_{i}+\phi\right)\cos\left(2\pi\omega\ T\right)+\cos\left(2\pi\omega\ x_{i}+\phi\right)\sin\left(2\pi\omega\ T\right)\right)=-\frac{A}{2\pi\omega}\cos\left(2\pi\omega\ T\right)\cos\left(2\pi\omega\ x_{i}+\phi\right)+\frac{A}{2\pi\omega}\sin{\left(2\pi\omega\ T\right)\sin\left(2\pi\omega\ x_{i}+\phi\right)}+\frac{A}{2\pi\omega}\cos\left(2\pi\omega\ x_{i}+\phi\right)-A\frac{T}{2}\sin\left(2\pi\omega\ x_{i}+\phi\right)-A\frac{T}{2}\cos\left(2\pi\omega\ T\right)\sin\left(2\pi\omega\ x_{i}+\phi\right)-A\frac{T}{2}\sin\left(2\pi\omega\ T\right)\cos\left(2\pi\omega\ x_{i}+\phi\right)=\left(-\frac{A}{2\pi\omega}\cos\left(2\pi\omega\ T\right)+\frac{A}{2\pi\omega}-A\frac{T}{2}\sin\left(2\pi\omega\ T\right)\right)\cos\left(2\pi\omega\ x_{i}+\phi\right)+\left(\frac{A}{2\pi\omega}\sin\left(2\pi\omega\ T\right)-A\frac{T}{2}+A\frac{T}{2}\cos\left(2\pi\omega\ T\right)\right)\sin\left(2\pi\omega\ x_{i}+\phi\right)\)(A8).
Inserting (A8) to the numerator of (A1), it becomes
\(\sigma_{L1}\left([x_{i},x_{i+1}]\right)=\left(-\frac{A}{2\pi\omega T}\cos\left(2\pi\omega\ T\right)+\frac{A}{2\pi\omega T}-\frac{A}{2}\sin\left(2\pi\omega\ T\right)\right)\cos\left(2\pi\omega\ x_{i}+\phi\right)+\left(\frac{A}{2\pi\omega T}\sin\left(2\pi\omega\ T\right)-\frac{A}{2}+\frac{A}{2}\cos\left(2\pi\omega\ T\right)\right)\sin\left(2\pi\omega\ x_{i}+\phi\right)\)(A9),
Now let's see the\(L1\left(\left|\varepsilon\right|\right)=L1\left(-\varepsilon\right)\)case. In this integrand, similarly
to (A3) becomes
\(\int_{x_{i}}^{x_{i+1}}{\left(\frac{\left(x_{i+1}-x\right)Asin(2\pi\omega\ x_{i}+\phi)+\left(x-x_{i}\right)Asin(2\pi\omega\ x_{i+1}+\phi)}{T}-Asin(2\pi\omega\ x+\phi)\right)\text{dx}}\)(A10),
so only the sign of the two main terms is exchanged. The derivation
using the steps through (A3) to (A9) are identical apart from the
consequent change of the sign, resulting in
\(\sigma_{L1}\left([x_{i},x_{i+1}]\right)=\left(\frac{A}{2\pi\omega T}\cos\left(2\pi\omega\ T\right)-\frac{A}{2\pi\omega T}+\frac{A}{2}\sin\left(2\pi\omega\ T\right)\right)\cos\left(2\pi\omega\ x_{i}+\phi\right)+\left(-\frac{A}{2\pi\omega T}\sin\left(2\pi\omega\ T\right)+\frac{A}{2}-\frac{A}{2}\cos\left(2\pi\omega\ T\right)\right)\sin\left(2\pi\omega\ x_{i}+\phi\right)\)(A11),
which is the same as (A9) apart from the opposite sign of the
coefficients of the sine and cosine functions of\(2\pi\omega\ x_{i}+\phi\).
Equations (A9) and (A11) provides a solution in a form of
\(\sigma_{L1}\left([x_{i},x_{i+1}]\right)=C\bullet\cos\left(2\pi\omega x_{i}+\phi\right)+S\bullet sin(2\pi\omega x_{i}+\phi)\)(A12).
By generalizing the solutions (A9) and (A11) in the form of (8) instead
of (A12), they are equivalent to equation (8) to (10) of the article.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig1/Fig1}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig2/Fig2}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig3/Fig3}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig4/Fig4}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig5/Fig5}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig6/Fig6}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig7/Fig7}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig8/Fig8}
\end{center}
\end{figure}
\selectlanguage{english}
\FloatBarrier
\end{document}