J. Appl. Geodesy 2022; aop
Sattar Isawi*, Harald Schuh, Benjamin Männel, and Pierre Sakic
Stability analysis of the Iraqi GNSS stations
https://doi.org/10.1515/jag-2022-0001
Received January 9, 2022; accepted February 20, 2022
Abstract: The Iraqi GNSS network was installed in 2005
with help from the USA and UK. The network consists of
seven GNSS stations distributed across Iraq. The network
GNSS data have been comprehensively analyzed in this
study; this, in turn, allowed us to assess the impact of var-
ious geophysical phenomena (e. g., tectonic plate motion
and Earthquakes) on its positional accuracy, stability, and
validity over time. We processed daily GPS data, spanning
over more than five years. The Earth Parameter and Or-
bit System software (EPOS.P8), developed by the German
Geoscience Research Center (GFZ), was used for data pro-
cessing by adopting the Precise Point Positioning (PPP)
strategy. The stacked time series of stations coordinates
was analyzed after estimating all modeled parameters of
deterministic and stochastic parts using the least-squares
technique. The study confirmed a slight impact of the re-
cent M 7.3 Earthquake on the Iraqi GNSS stations and con-
cluded that the stations are stable over the study period
(2013 up to 2018) and that the GNSS stations represent the
movement of the Arabian plate.
Keywords: CORS, time series analysis, power-law noise,
IGRS
1 Introduction
The crucial part of national infrastructure is establishing
an accurate spatial reference system that can be realized
by a reliable geodetic network of permanent control sta-
tions. These stations are referenced to a well-defined coor-
*Corresponding author: Sattar Isawi, Technische Universität Berlin,
Institute of Geodesy and Geoinformation Science, Strasse des 17.
Juni 135, 10623 Berlin, Germany, e-mail:
sattar.mb.isaw[email protected]u-berlin.de
Harald Schuh, Technische Universität Berlin, Institute of Geodesy
and Geoinformation Science, Strasse des 17. Juni 135, 10623 Berlin,
Germany; and GFZ German Research Centre for Geosciences,
Department of Geodesy, Telegrafenberg, 14473 Potsdam, Germany,
e-mail: harald.schuh@gfz-potsdam.de, ORCID:
https://orcid.org/0000-0001-5443-0370
Benjamin Männel, Pierre Sakic, GFZ German Research Centre for
Geosciences, Department of Geodesy, Telegrafenberg, 14473
Potsdam, Germany, e-mails: benjamin.maennel@gfz-potsdam.de,
pierre.sakic@gfz-potsdam.de, ORCID:
https://orcid.org/0000-0003-2938-1356 (B. Männel),
https://orcid.org/0000-0003-1770-0532 (P. Sakic)
dinate system to be the backbone of the nation’s geospatial
system (e. g., geodatabase and mapping systems).
Iraq followed this approach after the war in 2003
because reconstruction demanded implementing major
engineering projects. Moreover, the implementation of
these projects needs a stable and reliable geodetic basis.
Therefore, in 2005, Iraq established a new geodetic refer-
ence system called the Iraqi Geospatial Reference System
(IGRS). Since the establishment of IGRS, all new Iraqi en-
gineering projects and mapping systems were based on
it, specifically Oil & Gas, Transportation, Electrics, and
Agricultural projects. IGRS realized as a coincidence with
ITRF2000 at epoch 1997.0 but anchored to the Arabian
plate by five permanent GNSS stations [32]. In principle,
IGRS represents a snapshot of the Arabian plate’s dynamic
on which Iraq lies. Furthermore, IGRS is a dynamic ref-
erence frame that continuously moves with the Arabian
plate motion [10]. Besides, Iraq locates at the boundary be-
tween two different tectonic plates, Arabian and Eurasian
plates (Fig. 1). This area is deformable and continuously
unstable, as shown by several Earthquakes, for exam-
ple M 7.3 Iraq – Iran border earthquake on November 12,
2017 [26]. Quantifying the influence of such geophysical
processes on the stability of the Iraqi permanent GNSS sta-
tions necessitates an in-depth and comprehensive study
because of its significant impact on the present and future
geodetic surveying in Iraq.
The remarkable development of space geodetic tech-
niques (e. g., Global Navigation Satellite Systems GNSS)
made it easy to obtain continuous 3D positions for points
on Earth’s surface for very long periods of up to many
years. The obtained positions are characterized by suffi-
cient accuracy to monitor the GNSS station’s antenna posi-
tion change over time. This change indirectly refers to var-
ious physical processes in the Earth system (e. g., tectonic
plate motion).
Based on this proposition, several researchers stud-
ied the movement of the Arabian plate by adopting space
geodetic techniques. For example, [3] studied the Arabian
plate motion by estimating the Riyadh SLR station’s posi-
tion and velocity. Laser-ranging observations of 14 years
from 20 SLR stations to LAGEOS-1 and LAGEOS-2 have
been analyzed in this study. They concluded that the Ara-
bian plate is countering a north-east motion with an an-
nual linear rate of 42.9 mm/year. [11] processed ten years
of daily GPS data from two Iraqi GNSS stations. ISER, lo-
cated in Erbil, and ISNA, located in Najaf. The GITSA anal-
Open Access. © 2022 Isawi et al., published by De Gruyter. This work is licensed under the Creative Commons Attribution 4.0 International
License.
2| S. Isawi et al., Stability analysis of the Iraqi GNSS stations
Figure 1: The spatial distribution of the Iraqi GNSS stations. The black dashed lines are the boundary of tectonic plates. The star symbol
represents the epicentre location of the M 7.3 Earthquake that hits Iraq – Iran border on November 12, 2017.
ysis software was used in this study to build up and ana-
lyze the time series of station coordinates. The study con-
cluded that the estimated velocity vectors for the ana-
lyzed stations were 38 mm/year and 40 mm/year, respec-
tively.
The estimation of linear rate (mainly driven by plate
motion) is the main focus of these studies. However, this
estimation needs an improvement to the adopted func-
tional model. The periodic effects (e. g., annual and semi-
annual deformation), sudden effects (e. g., offset due to
antenna change or co-seismic deformation), or even post-
seismic deformation have to be modeled and extracted
from the time series signal [4]. Although all common de-
terministic parts would be modeled for the time series, dif-
ferent noise sources are still present in its signal due to
satellite orbit error, atmospheric delays, and clock insta-
bility. On the other hand, the assumption that the only
white noise (uncorrelated noise) contributes to the noise
model will lead to underestimated uncertainty of linear ve-
locities by factors from three to six times [33], or even five
to 11 times [20].
The Iraqi GNSS stations will be extensively analyzed
in this study to assess the impact of various geophysical
processes (e. g., tectonic plate motion and Earthquakes) on
their accuracy, stability, and validity over the study period.
The determination of coordinate time series characteris-
tics (trend, annual signals, offsets, and spectra) will be fur-
ther described and discussed. Besides, this study will in-
vestigate the potential contribution of different noise types
(e. g., power-law noise) to the noise model. A Matlab code
comprises functions written to implement this study in
both the frequency and time domain.
2 Data
Five continuously operating GNSS reference stations have
been selected out of seven Iraqi GNSS stations for this
S. Isawi et al., Stability analysis of the Iraqi GNSS stations | 3
Figure 2: GPS antennas attached to the monuments. ISER, ISKU (top)
and ISBA, ISNA (middle), and ZAXO station (bottom). Source NGS at
https://www.ngs.noaa.gov/.
study. The selected GNSS stations are ISNA, ISKU, ISBA,
ISER, and ZAXO. ISBS and ISSD stations were excluded be-
cause both of them were non-operational during the inves-
tigated period. The spatial distribution of the Iraqi GNSS
stations is shown in Fig. (1). All selected GNSS stations
were equipped with dual-frequency GNSS receivers as well
as geodetic antenna types. However, all of them have non-
geodetic monuments because of the antennas mounted on
top of buildings (Fig. 2). Station ISBA is part of the Inter-
national GNSS Service (IGS) network [12], whereas the re-
maining stations are a part of the National Geodetic Survey
(NGS) global network.
The daily observation (raw data) of all GNSS stations
were downloaded from the NGS public archive at ftp://
geodesy.noaa.gov/cors/rinex/. The number of daily obser-
vation files for each station is shown in Fig. (3). The GNSS
station’s meta-data (log files) needed to track the station
history were also obtained from the NGS database system
(Table 1). The data included a comprehensive record of sta-
Figure 3: Number of daily observation files per station.
tion firmware and hardware changes (e. g., GNSS receiver
and antenna type). A pre knowledge in station history is
required in this study because the change in station hard-
ware and software could lead to offsets in the position time
series [27].
2.1 GPS daily solutions
The software package EPOS.P8, developed at GFZ, was
used to process the GPS data. Satellite orbits, clock cor-
rections, and Earth Orientation Parameters (EOP) were
provided by a GFZ internal reprocessing based on the
operational final solution [18]. The products were given
in the IGS14 reference frame. Our final solution of the
GPS daily coordinates for each station was estimated fol-
lowing the Precise-Point-Positioning (PPP) methodology
adopted in GFZ AC. The processing details (Table 11) and
the flow diagram (Fig. 8) of PPP processing are shown in
Appendix A.
3 Methods
Firstly, the 3-D Cartesian coordinates Xi,Yi, and Ziob-
tained from GNSS data processing were transformed into
topocentric North, East, and Up. The coordinate trans-
formation from the right-hand system Earth Centered
Table 1: Log files information of selected GNSS Stations.
ISNA ISKU ISBA ISER ZAXO
City Najaf Kut Baghdad Erbil Dohuk
Sat. System GPS+GLO GPS+GLO GPS+GLO GPS+GLO GPS
GPS Receiver TRM NETR9 TRM TRM NETR5 TRM NETR5 TPS NET-G3A
Ant. Type TRM57971 TRM57971 TRM41249 TRM57971 TPSCR.G5
Monum. Desc. Tribr. adap. Trib. adap. Trib. adap. Trib. adap. Trib. adap.
4| S. Isawi et al., Stability analysis of the Iraqi GNSS stations
Earth Fixed (ECEF) into the topocentric left-hand system
has been implemented according to [24]. The variance-
covariance matrix SXYZ of GNSS daily solutions have also
been transformed into the topocentric system by applying
a variance-covariance propagation law [25].
3.1 Functional model
The time series for individual topocentric components
(North, East, and Up) was built up based on the GNSS daily
solutions ytiat uniformly spaced time ti(i=1,...,n, where n
is the number of daily solutions). The stacked time series
have been modeled (approximated) as the sum of deter-
ministic parts, by considering the initial site position x0,
and linear rate vx(interpreted as the GNSS station veloc-
ity), as well as the periodic signals of annual and semi-
annual motion. The functional model reads:
yti=x0+vxti+asin(2πti)+bcos(2πti)+csin(4πti)
+dcos(4πti)+εi(1)
with aand bare the coefficients of the annual periodic
motion, while cand dare the coefficients of semi annual
motion. When there are in addition moffsets (jumps) in
the coordinate time series, the functional model can be ex-
panded to [6]:
yti=x0+vxti+asin(2πti)+bcos(2πti)+csin(4πti)
+dcos(4πti)+m
∑
j=1
fjh(ti−tfj)+εi(2)
where fjis the offset’s magnitude, tfj is the offset’s epoch
with hdenoting the Heaviside function, and εis the model
error. With apriori known offset epochs tfj, the linear model
of observation equations in Eq. (2) is solved using the least-
squares method:
x=(ATSyy−1A)−1ATSyy−1y,(3)
where the vector of adjusted unknowns
xconsists of the
following parameters:
x=[
x0
vx
a
b
c
d
f]T,(4)
Syy is the variance-covariance matrix of the observation
vector y,
y=[yt1,...,ytn]T,(5)
and, the design matrix Ais given as:
A=[[[[
1t1sin(2πt1)cos(2πt1)sin(4πt1)cos(4πt1)h(t1−tfn)
.
.
..
.
..
.
..
.
..
.
..
.
..
.
.
1tnsin(2πtn)cos(2πtn)sin(4πtn)cos(4πtn)h(tn−tfm)]]]]
(6)
Table 2: Summary of coordinate data.
Site Data Start epoch End epoch Outliers
# # %
ISBA
N 1743 2013.00 2018.99 134 7.6
E 1743 2013.00 2018.99 116 6.6
U 1743 2013.00 2018.99 145 8.3
ISNA
N 1327 2013.00 2018.99 89 6.7
E 1327 2013.00 2018.99 88 6.6
U 1327 2013.00 2018.99 51 3.8
ISKU
N 1591 2013.00 2018.99 70 4.3
E 1591 2013.00 2018.99 78 4.9
U 1591 2013.00 2018.99 100 6.2
ISER
N 1837 2013.00 2018.99 284 15.4
E 1837 2013.00 2018.99 267 14.5
U 1837 2013.00 2018.99 155 8.4
ZAXO
N 1090 2015.00 2018.99 110 10
E 1090 2015.00 2018.99 69 6.3
U 1090 2015.00 2018.99 65 5.9
3.2 Outlier detection and removal
Outliers have been detected and removed from each co-
ordinate time series during the preliminary model fitting.
The median and interquartile range approach [21] were im-
plemented in this study. The median and IQR statistics are
calculated for the residuals on a sliding window of one
year centered on each time series element to consider the
gradual decrease in the data variance associated with the
increase in the time series length. Outliers were defined as
having:
|vi−median|>1.5×IQR (7)
where IQR defines the middle 50 % of the sample data, and
viis the ith least square residual.
Summary of removed outliers for each station time se-
ries is presented in Table (2).
4 Noise characteristics and
analysis
As with many other geophysical phenomena, noise in co-
ordinate time series can be described as a power law pro-
cess of the form [2, 16]:
S. Isawi et al., Stability analysis of the Iraqi GNSS stations | 5
Pv(f)=P0(f
f0)k
(8)
where fis the temporal frequency, P0and f0are normal-
izing constants, and kis the spectral index that defines a
slope of the best fit line to the spectra presented in log-log
space [33]. A smaller spectral index implies a more corre-
lated process and a more relative power at lower frequen-
cies. Therefore such a natural process is characterized by
negative indices [13]. According to [17], the spectral index
of different geophysical phenomena is classified into two
groups. The first group is called fractional Brownian mo-
tion with spectral index −3<k< −1. The noise processes
in this interval are non-stationary including random walk
with k= −2 or(Px∝1/f2), while the second group termed
fractional Gaussian with spectral index −1<k< +1, where
the noise processes in this interval are stationary (inde-
pendent of time) including the special case of uncorrelated
white noise with k=0 or (Px∝1/f0). The special case
k= −1 or (Px∝1/f1) is known as “flicker noise”. Plenty
of studies showed that the spectral index has particular
importance because it is a good indicator for characteriz-
ing the noise source [14, 22]. Based on the preceding, it is
necessary to study the post-fit residuals in-depth to decide
which colored (correlated) noise type has to be considered
in the noise model.
Two approaches are described in the literature to clas-
sify and quantify noise in the coordinate time series. The
first approach’s implementation relies on time series anal-
ysis in the frequency domain (spectral analysis). In con-
trast, the second one relies on the time series analysis in
the time domain. In this study, we restricted ourselves to
the second approach, but in addition, only estimating the
spectral index in the frequency domain.
4.1 Noise analysis
The time series signal and its noise do not behave in the
same way when transformed and plotted as a power spec-
trum. This property can be of great benefit for noise anal-
ysis in the frequency domain. Based on this concept, the
time series of post-fit residuals have been transformed into
the frequency domain for further analysis.
Unfortunately, equipment malfunction, interruptions
in the communications network, or even outlier removal
lead to missing data in the coordinate time series. There-
fore, a Lomb-Scargle periodogram introduced by [15] and
used by [33] has been applied to solve this issue.
According to [33], the time series of post-fit residuals
must be long enough and comprise data with a range of fre-
quencies that make the power spectra well approximated.
Hence, from the practical aspect, the spectral index of the
power-law process can be estimated by fitting a straight
line to these frequencies to avoid biasing the regression es-
timate with the predominance of white noise at high fre-
quencies and decreasing the impact of outliers at lower
frequencies. The functional model of power-law process in
Eq. (8) can be reformulated as:
Pvi(f)=P0
f0kfi
k(9)
substitution leads to
W0=P0
f0k(10)
Taking the log of both sides yields
log(Pvi(f)) =log(W0)+klog(fi)(11)
where all parameters are defined as in Eq. (8). To estimate
the slope of the fitted line (spectral index k) to the spec-
tra at log-log space, the least-squares method was used.
Therefore, Eq. (3) has been applied with the design matrix
of the form:
A=[[[[[[[[[[
1 log(f1)
1 log(f2)
. .
. .
. .
1 log(fN)
]]]]]]]]]]
(12)
The spectral index has been estimated for each time se-
ries. The proper noise model was founded on a combina-
tion of white and flicker noises. Since the amplitudes of
both noises have to be known in the adopted noise model,
the contribution for each noise has been estimated in this
study using the Maximum likelihood method presented
in [30]. In this approach, the amplitudes (variances) were
estimated by explicitly maximizing the likelihood function
for the total variance σ2
σ2=
vTC−1
v
N(13)
where C is the covariance matrix of post-fit residuals, N is
the number of daily solutions, and
v is the vector of post-
fit residuals. According to [30], the covariance matrix C can
be decomposed into two components with amplitudes σ2
w
and σ2
kbeing the variance of white and power-law noises
respectively such that:
C=σ2
wIw+σ2
kJk(14)
6| S. Isawi et al., Stability analysis of the Iraqi GNSS stations
where Iwis the identity matrix, and Jkis the covariancema-
trix of power-law noise, and both of them have the size of
N×N. The elements of matrix Jkwith k= −1 (flicker noise)
have been approximated in this study according to [33]:
J-1 ≡{
{
{
{
{
{
{
{
{
{
{
{
{
{
{
{
{
(3
4)2×2ti=tj
(3
4)2×(2−
log(ti−tj)
log(2)+2
12 )ti= tj
(15)
The approximated covariance matrix in Eq. (15) is not ex-
actly the same as that derived by the method of fractional
difference, introduced and used by [9], [28], and [30]. The
only difference is in the scaling [29].
5 Results and discussion
High precision coordinate time series have been produced
after processing GNSS data and cleaning for outliers. The
most common deterministic and stochastic parts were
modeled and described in Eq. (2). Various model param-
eters were estimated and assessed. In this section, we
provide further details and discuss all estimated param-
eters.
5.1 Stochastic part
It is necessary to investigate the post-fit residuals further to
build up a proper noise model. Hence, real uncertainties
can be assigned to the estimated long-term rate (velocity)
and other parameters in the underlying functional model.
A typical example of estimated power spectral density and
its associated spectral index is given for the station ISKU.
The station observations span over more than five years of
GNSS data as presented in Fig. (4). At this figure, the noise
power is almost constant for up to cross-over frequency
at (1/12 day) which implies the predominant white noise
at high frequencies. In contrast, the power increases from
a cross-over frequency to the left, indicating the predomi-
nant time-dependent noise at low frequencies.
Furthermore, due to errors in the solid Earth tide
model, or in the indirect tidal errors in the satellite or-
bits [1], a noise peak at the frequency (1/13.6 day) is ob-
served in this site signal, and at most analyzed time series
as well. However, the estimation of spectral indices was
done for frequencies more minor than (1/14 day) to prevent
biasing the estimation with this peak and exclude the pre-
dominant white noise.
Figure 4: The Lomb-Scargle periodogram for the station ISKU. The
solid lines (magenta) are the best fit in log-log space from frequen-
cies less than the cross-over frequency 1/(14 days).
Figure 5: Estimated spectral indices.
The spectral analysis shows that the estimated spec-
tral indices (Fig. 5) range from −0.20 ±0.05 to −0.50 ±0.04
for the north direction and −0.48 ±0.04 to −0.59 ±0.04
for the east direction, while the indices for up direction
range from −0.37 ±0.08 to −0.57 ±0.03. Furthermore,
the estimated spectral indices for all time series are rang-
S. Isawi et al., Stability analysis of the Iraqi GNSS stations | 7
Table 3: Estimated white + flicker noise amplitudes.
Site Spectral Index W.Noise
±σ(mm)
% F.Noise
±σ(mm)
%
ISBA
N−0.46 ±0.06 0.41 ±0.03 5 3.16 ±0.23 95
E−0.79 ±0.06 0.65 ±0.16 6 4.32 ±0.02 94
U−0.68 ±0.05 1.45 ±0.02 8 8.45 ±0.08 92
ISNA
N−0.13 ±0.06 0.19 ±0.02 1 3.31 ±0.29 99
E−0.46 ±0.05 0.31 ±0.24 2 4.09 ±0.02 98
U−0.47 ±0.06 1.14 ±0.01 5 8.69 ±0.10 95
ISKU
N−0.41 ±0.06 0.58 ±0.03 6 3.85 ±0.19 94
E−0.59 ±0.06 0.81 ±0.15 8 4.85 ±0.03 92
U−0.54 ±0.06 1.29 ±0.01 6 9.21 ±0.08 94
ISER
N−0.51 ±0.05 0.76 ±0.03 10 4.01 ±0.17 90
E−0.73 ±0.05 0.69 ±0.15 6 4.63 ±0.02 94
U−0.54 ±0.05 1.03 ±0.01 4 9.10 ±0.08 96
ZAXO
N−0.24 ±0.05 0.25 ±0.02 1 3.88 ±0.32 99
E−0.51 ±0.06 0.58 ±0.25 5 4.65 ±0.03 95
U−0.33 ±0.06 1.36 ±0.02 6 9.22 ±0.12 94
ing within fractional Gaussian part (−1<k<0). This
means the power-law noise being the main contributor to
the character of stochastic part, confirming the findings of
[33], [20], and [30]. Hence, the noise model of each time
series was determined accordingly.
Further to preceding spectral analysis, the adopted
noise model identified previously was analyzed. This anal-
ysis included quantifying the contribution of white and
flicker noises by estimating the amplitudes (variances) in
the time domain. The analysis shows a predominant con-
tribution for flicker noise to the total noise model. This con-
tribution ranges from 90 % to 99 %, while the white noise
contribution ranges from 1 % to 10 %. The initial findings
suggested to perform the data fitting once more for all co-
ordinate time series, but in this case, with a presumed
stochastic model mainly based on the estimated contribu-
tion for each noise type. The estimated amplitudes are pro-
vided in Table (3). This Table shows that the estimated am-
plitudes of flicker noise are very consistent for all coordi-
nate time series. The same also applies to white noise, indi-
cating that a similar noise type contaminates all analyzed
coordinate time series.
5.2 Regression analysis
The completion of the noise analysis created a solid
ground to implement the regression analysis. This, in turn,
Figure 6: Estimated offsets due to M 7.3 Earthquake.
allowed us to adopt a realistic stochastic model to pro-
ceed with the weighted least-squares estimation for the
data model in Eq. (2). Various model parameters were es-
timated and assessed in this analysis, for example, the
time series jump (mainly driven by Earthquakes), trend
(plate motion), and the annual and semi-annual variation
of coordinate time series. In this section, we discuss in de-
tail the outcomes of the regression analysis we have per-
formed.
5.2.1 Discontinuities due to Earthquake
We evaluated the impact of the recent M 7.3 Earthquake on
the stability of the network stations. The evaluation has
been done simultaneously with testing the slight potential
effect of receiver hardware and firmware change on the co-
ordinate time series continuity. The findings of this evalua-
tion are presented in Fig. (6). This figure shows that all sta-
tions displaced in the up direction except the station ISER
which is only displaced in the north direction. ZAXO sta-
tion was relatively more influenced (4.5 mm offset at up di-
rection) than the remaining stations whose offsets ranged
from 2.4 mm to 3.40 mm at up direction. Details of all ob-
served offsets caused by Earthquake and firmware change
are provided in Table (4).
5.2.2 Station velocities
The outcomes of the regression analysis show very similar
horizontal velocities for all stations (Table 5). It is an ex-
pected result as all stations are close to each other within
the rigid part of the Arabian plate. The estimated veloci-
ties after removing all observed offsets (jumps) range from
25.0±0.2 mm/year to 27.7±0.1 mm/year for the north
direction and 22.1±0.2 mm/year to 24.9±0.1 mm/year
for the east direction. The velocity vectors depicted in
8| S. Isawi et al., Stability analysis of the Iraqi GNSS stations
Table 4: Observed Offsets.
Site Axis Offset (mm) Epoch Caused by
ISBA E 1.70 ±0.25 2017.8657 M 7.3 earthquake
ISBA U 2.50 ±0.47 2017.8657 M 7.3 earthquake
ISBA U −3.50 ±0.65 2015.2054 firmware change
ISNA E −1.01 ±0.43 2017.7534 firmware change
ISNA U 2.40 ±0.95 2017.8657 M 7.3 earthquake
ISKU N −1.20 ±0.30 2017.8657 M 7.3 earthquake
ISKU U 3.40 ±0.71 2017.8657 M 7.3 earthquake
ISKU U −1.23 ±0.60 2015.1643 firmware change
ISKU U 1.97 ±0.68 2016.9890 firmware change
ISER N −1.10 ±0.25 2017.8657 M 7.3 earthquake
ISER N 2.40 ±0.23 2014.0958 unknown
ISER N 1.30 ±0.21 2014.6876 firmware change
ISER E 1.40 ±0.27 2014.0958 unknown
ISER E −1.20 ±0.55 2014.6876 firmware change
ZAXO E 1.30 ±0.35 2016.1041 firmware change
ZAXO U −1.70 ±0.69 2016.1041 firmware change
ZAXO U 4.50 ±0.79 2017.8657 M 7.3 earthquake
Table 5: Estimated staions velocities.
Site Velo. (mm/yr) σ(mm/yr) Velo. (mm/yr) σ(mm/yr)
White+Flicker White-only
ISBA
N 27.67 ±0.13 27.79 ±0.03
E 24.11 ±0.16 24.07 ±0.04
U−0.92 ±0.24 −1.25 ±0.08
ISNA
N 27.27 ±0.10 27.20 ±0.04
E 24.85 ±0.12 24.73 ±0.05
U−1.98 ±0.23 −1.93 ±0.09
ISKU
N 26.18 ±0.15 26.25 ±0.04
E 23.26 ±0.18 23.29 ±0.05
U−3.30 ±0.23 −3.51 ±0.09
ISER
N 25.03 ±0.18 24.96 ±0.04
E 22.14 ±0.17 21.87 ±0.05
U−4.92 ±0.20 −5.05 ±0.08
ZAXO
N 26.27 ±0.14 26.82 ±0.08
E 22.82 ±0.21 23.00 ±0.09
U−1.39 ±0.32 −2.02 ±0.18
Fig. (7) shows that all stations experienced a north-east
movement, which mainly represents a snapshot of the
Arabian plate motion in Iraq. In contrast, our study
found that the vertical velocities −4.9±0.2 mm/year and
−3.3±0.2 mm/year for ISER and ISKU respectively (see
Appendix B and C), are remarkable vertical velocities and
worthy of attention and assessment. The most realistic ex-
planation to the observed velocities might be the effect
Figure 7: Locational map of the M 7.3 Earthquake epicenter (gray
star) and the depicted velocity vectors (black arrows) of the Iraqi
GNSS stations.
driven by non-tectonic signals. This affection can cause
deformation, for example, groundwater withdrawal pro-
duced long-term subsidence [5]. The less likely but inter-
esting explanation is that all network stations have a non-
geodetic monument because the antennas were mounted
on buildings (Fig. 2). Consequently, the buildings them-
selves could be subject to deformation. On the other hand,
we estimated the velocity rates for all stations two times
by applying two different stochastic models to emphasize
the importance of the implemented stochastic model in
this study. First, the velocity rate was estimated using only
a white noise-based stochastic model, while at the sec-
ond time, a stochastic model comprised white noise plus
power-law noise (flicker noise). Table (5) shows the esti-
mated velocities and their uncertainties in both cases. As
expected, the resulting uncertainties for case one (white
noise-based stochastic model) are much smaller than the
second case. The underestimation by almost one order of
magnitude confirms the findings of [33], and [21].
5.2.3 Annual variations
Annual and semi-annual variations have been observed
in all coordinate time series, at both horizontal and verti-
cal components. The annual variations range from 0.5±
0.1 mm to 2.5±0.1 mm for the north direction and 0.2±
0.1 mm to 3.9±0.1 mm for the east direction, while in the
up direction the situation is quite different. The annual
variations at vertical components range from 7.8±0.2 mm
to 12.8±0.1 mm. Although these results are in good agree-
S. Isawi et al., Stability analysis of the Iraqi GNSS stations | 9
Table 6: Estimated annual signals.
Site Amplitude
(mm)
Uncertainty
(mm)
Phase
(deg)
Uncertainty
(deg)
ISBA
N 1.49 ±0.07 −118.32 ±0.002
E 2.02 ±0.09 −080.95 ±0.002
U 12.85 ±0.05 −112.74 ±0.006
ISNA
N 0.88 ±0.08 −163.57 ±0.006
E 0.24 ±0.10 021.99 ±0.007
U 10.98 ±0.20 −112.11 ±0.008
ISKU
N 2.46 ±0.08 082.00 ±0.001
E 1.00 ±0.11 −089.51 ±0.001
U 10.14 ±0.20 −117.66 ±0.008
ISER
N 0.50 ±0.08 −141.63 ±0.006
E 2.61 ±0.10 −085.92 ±0.001
U 7.72 ±0.18 −108.17 ±0.006
ZAXO
N 0.73 ±0.10 −069.68 ±0.004
E 3.86 ±0.11 104.46 ±0.003
U 12.39 ±0.24 −134.17 ±0.015
ment with each other, the loading corrections, as sug-
gested for example by [19] were not applied in process-
ing GNSS data. Thus, the height variations are much larger
than the variations in the horizontal components. Further-
more, the local hydrological loading causes height varia-
tions of (10–30 mm) with a predominant annual period [8].
For example, the Tigris and Euphrates rivers and the sev-
eral dams distributed at middle and north of Iraq. These
rivers and dams are usually full in the winter and are rel-
atively empty in the summer season. Besides, the thermal
expansion of the station monuments could contribute to
this oscillation [23]. The stations monuments are exposed
to a significant temperature variation due to the temper-
ature change in Iraqi’s summer and winter season. The
average daily temperature in winter reaches 6 ∘C, while
48 ∘C in summer [31]. Hence, the station monument may
experience a thermal expansion because the GNSS anten-
nas are mounted on galvanized iron rods at top of the
buildings (Fig. 2). Details of estimated annual and semi-
annual variations are presented in Tables (6) and (7), re-
spectively. The post-fit residuals obtained after removing
the annual variations from all coordinate time series are
presented in Table (8). This table shows very similar resid-
uals for all analyzed time series, indicating a similar noise
type.
Table 7: Estimated semi-annual signals.
Site Amplitude
(mm)
Uncertainty
(mm)
Phase
(deg)
Uncertainty
(deg)
ISBA
N 1.04 ±0.05 141.76 ±0.002
E 1.19 ±0.05 130.54 ±0.001
U 0.95 ±0.01 117.74 ±0.003
ISNA
N 0.81 ±0.07 130.18 ±0.015
E 0.81 ±0.33 107.73 ±0.006
U 0.74 ±0.01 047.50 ±0.005
ISKU
N 1.17 ±0.04 172.98 ±0.001
E 0.74 ±0.08 155.97 ±0.001
U 1.40 ±0.03 061.36 ±0.004
ISER
N 0.58 ±0.09 178.21 ±0.013
E 0.80 ±0.03 151.74 ±0.001
U 1.98 ±0.04 085.68 ±0.001
ZAXO
N 0.75 ±0.11 113.82 ±0.001
E 0.80 ±0.02 155.34 ±0.003
U 2.40 ±0.04 048.98 ±0.014
Table 8: Post-fit residuals scatter.
Site N (mm) E (mm) U (mm)
ISBA 1.90 2.60 5.10
ISNA 1.90 2.40 5.10
ISKU 2.30 2.90 5.40
ISER 2.40 2.70 5.30
ZAXO 2.20 2.70 5.50
6 Validating the results
The estimated velocities of the GNSS stations in this study
were validated with two independent studies. However,
the two studies did not account for the discontinuities (off-
sets) and the annual variation of the coordinate time se-
ries. The velocity vectors estimated by [3] using Satellite
Laser Ranging for the SLR-station at Kingdom of Saudi Ara-
bia, also located at the stable part of the Arabian plate,
agrees very well with the velocity rates observed in this
study (Table 9). The second validation (Table 10) shows,
that the velocities of the current study are 4 mm/year less
than those estimated by [11] that was obtained from pro-
cessing ten years of daily observation data for two GNSS
stations in mid and north of Iraq.
10 | S. Isawi et al., Stability analysis of the Iraqi GNSS stations
Table 9: Summary of first comparsion (GNSS vs. SLR).
Study Site Estimated Velocity Velocity. Vec. (mm/yr) Azimuth Technique
N E U
⋆Riyadh 29.10 31.60 1.90 42.90 N47°24´E SLR
⋆⋆ ISBA 27.67 24.11 −0.92 36.70 N41°03´E GNSS
ISNA 27.27 24.85 −1.98 36.70 N42°20´E GNSS
ISKU 26.18 23.26 −3.30 35.02 N41°36´E GNSS
ISER 25.03 22.14 −4.92 33.41 N41°29´E GNSS
ZAXO 26.27 22.82 −1.39 34.79 N40°58´E GNSS
⋆Alothman and Schillak [3] ⋆⋆This study
Table 10: Summary of second comparsion.
Study Site Velocity Vector
(mm/yr)
Azimuth Space Technique
⋆ISER 38.00 N40°48´E GNSS
ISNA 40.00 N45°06´E GNSS
⋆⋆ ISER 33.41 N41°29´E GNSS
ISNA 36.70 N42°20´E GNSS
⋆Jasim et al. [11] ⋆⋆ This study
7 Conclusion
We assessed the stability of Iraq’s core GNSS stations
based on the processing and interpretation of the GNSS
daily data from 2013 up to 2018, derived for five perma-
nent GNSS stations across Iraq. Determining coordinate
time series characteristics (trend, annual signals, discon-
tinuities, and spectra) were comprehensively discussed.
The spectral analysis showed that the estimated spec-
tral indices for all time series fall within the fractional
Gaussian part (−1<k<0). This, in turn, refers to the
power-law noise being the main contributor to the stochas-
tic part. Besides, the estimated annual signals with ampli-
tudes of around 10 mm were well observed in the up com-
ponent for all time series. The signals agreed with each
other and reflected the impact of different geophysical pro-
cesses in the Earth system. The study concluded, that the
source of such annual signals mainly refers to the near-
surface hydrologic loading and the seasonal temperature
variation.
The study shows a slight impact of the recent M 7.3
Earthquake on the Iraqi GNSS stations. The maximum off-
set does not exceed five milimeters. Moreover, our anal-
ysis indicates that several offsets caused by equipment
firmware change occurred; and were removed in the coor-
dinate time series.
The determination of station velocities included a
comparison between different noise assumptions. The
study reveals the benefit of choosing a suited noise model
(e. g., the combination of white and power-law noises)
on the derived velocities. The estimated stations veloci-
ties range from 25 ±0.1 mm/year to 27.6±0.1 mm/year
for north direction and 22.1±0.1 mm/year to 24.8±
0.1 mm/year for east direction. Subsidence rates of
−4.9 mm/year and −3.3 mm/year are observed for stations
ISER and ISKU, respectively.
A comparison was made with two former studies to
validate the results of this study. The validation showed
very good agreement with the results of [3, 11], and the
study confirmed that the Iraqi GNSS stations are stable
over the study period and that the stations represent the
Arabian plate’s motion.
S. Isawi et al., Stability analysis of the Iraqi GNSS stations | 11
Appendix A. PPP processing informations
Table 11: Basic models and parameters of PPP processing.
Modeling and a priori information
Observations ionosphere free linear combination formed by undifferenced GPS observation, sampled with 300 s,
elevation cut-off 5°
A priori products GFZ SP3 orbits, clock corrections, Earth orientation parameters [18]
Tropospheric correction empirical model GPT2 mapped with GMF [7]
Ionospheric correction first-order effect considered with ionosphere-free linear combination
GNSS phase center igs14.atx
Estimated parameters
Coordinates estimated daily in IGS14 reference frame introduced via the orbit products
Troposphere 25 zenith delays per day; two gradient pairs per station and day
Receiver clock pre-eleminated every epoch
GNSS ambiguities estimated, without ambiguity fixing
Figure 8: Flow diagram of PPP processing via EPOS.P8 software.
12 | S. Isawi et al., Stability analysis of the Iraqi GNSS stations
Appendix B. Position time series of site ISKU
Figure 9: The Top, middle, and bottom are the position time series for north, east, and up respectively. The daily solutions (black dots) of the
site ISKU are modeled (red lines) with a linear trend. The blue lines refer to the removed related offsets.
S. Isawi et al., Stability analysis of the Iraqi GNSS stations | 13
Appendix C. Position time series of site ISER
Figure 10: The Top, middle, and bottom are the position time series for north, east, and up respectively. The daily solutions (black dots) of
the site ISER are modeled (red lines) with a linear trend. The blue lines refer to the removed related offsets.
14 | S. Isawi et al., Stability analysis of the Iraqi GNSS stations
References
[1] Abraha, K.E., F.N. Teferle, A. Hunegnaw, and R. Dach. 2017.
GNSS related periodic signals in coordinate time-series from
precise point positioning. Geophysical Journal International
208(3): 1449–1464. DOI:10.1093/gji/ggw467.
[2] Agnew, D.C. 1992. The time-domain behavior of power-law
noises. Geophysical research letters 19(4): 333–336.
[3] Alothman, A. and S. Schillak. 2014. Recent results for the
arabian plate motion using satellite laser ranging observations
of riyadh SLR station to LAGEOS-1 and LAGEOS-2 satellites.
Arabian Journal For Science and Engineering 39(1): 217–226.
DOI:10.1007/s13369-013-0823-7.
[4] Altamimi, Z., P. Rebischung, L. Métivier, and X. Collilieux.
2016. ITRF2014: A new release of the International Terrestrial
Reference Frame modeling nonlinear station motions. Journal
of Geophysical Research: Solid Earth 121(8): 6109–6131.
DOI:10.1002/2016JB013098.
[5] Bawden, G.W., W. Thatcher, R.S. Stein, K.W. Hudnut, and
G. Peltzer. 2001. Tectonic contraction across los angeles after
removal of groundwater pumping effects. Nature 412(6849):
812–815. DOI:10.1038/35090558.
[6] Bevis, M. and A. Brown. 2014. Trajectory models and reference
frames for crustal motion geodesy. Journal of Geodesy 88(3):
283–311. DOI:10.1007/s00190-013-0685-5.
[7] Böhm, J., B. Werl, and H. Schuh. 2006. Troposphere mapping
functions for gps and vlbi from ecmwf operational analysis
data. Journal of geophysical research 111(B2): B02406.
[8] Dill, R. and H. Dobslaw. 2013. Numerical simulations
of global-scale high-resolution hydrological crustal
deformations. Journal of Geophysical Research: Solid Earth
118(9): 5008–5017. DOI:10.1002/jgrb.50353.
[9] Hosking, J. 1981. Fractional differencing. biometrika 68:
165–176. Mathematical Reviews (MathSciNet): MR614953
Zentralblatt MATH 464.
[10] International Association of Oil & Gas Producers. 2015, may.
Geodetic referencing in Iraq. Technical report.
[11] Jasim, Z.N., O.Y. Alhamadani, and M.U. Mohammed. 2018.
Investigation the plate motion using operating refere.
International Journal of Civil Engineering 9(13).
[12] Johnston, G., A. Riddell, and G. Hausler. 2017. The
international GNSS service, In Springer handbook of global
navigation satellite systems, 967–982. Springer.
[13] Kenyeres, A. and C. Bruyninx. 2009. Noise and periodic terms
in the EPN time series, In Geodetic Reference Frames, 143–148.
Springer. DOI:10.1007/978-3-642-00860-3_22.
[14] Langbein, J., F. Wyatt, H. Johnson, D. Haman, and P. Zimmer.
1995. Improved stability of a deeply anchored geodetic
monument for deformation monitoring. Geophysical Research
Letters 22(24): 3533–3536. DOI:10.1029/95GL03325.
[15] Lomb, N. and J. Scargle. 1976. Further development and
properties of the spectral analysis by least-squares.
Astrophysics and Space Science 39: 447.
[16] Mandelbrot, B.B. 1983. The fractal geometry of nature/revised
and enlarged edition. New York, WH Freeman and Co., 1983,
495 p.
[17] Mandelbrot, B.B. and J.W. Van Ness. 1968. Fractional brownian
motions, fractional noises and applications. SIAM review
10(4): 422–437. DOI:10.1137/1010093.
[18] Männel, B., A. Brandt, T. Nischan, A. Brack, P. Sakic, and
M. Bradke. 2020. GFZ final product series for the International
GNSS Service (IGS). DOI:10.5880/GFZ.1.1.2020.002.
[19] Männel, B., H. Dobslaw, R. Dill, S. Glaser, K. Balidakis,
M. Thomas, and H. Schuh. 2019. Correcting surface loading
at the observation level: impact on global gnss and vlbi
station networks. Journal of Geodesy 93(10): 2003–2017.
DOI:10.1007/s00190-019-01298-y.
[20] Mao, A., C.G. Harrison, and T.H. Dixon. 1999. Noise in GPS
coordinate time series. Journal of Geophysical Research: Solid
Earth 104(B2): 2797–2816. DOI:10.1029/1998JB900033.
[21] Nikolaidis, R. 2002, January. Observation of geodetic and
seismic deformation with the Global Positioning System. Ph. D.
thesis, University of California, San Diego.
[22] Nistor, S. and A. Buda. 2014. Time series analysis in the
presence of coloured noise. Journal of Applied Engineering
Science 4(17): 2.
[23] Romagnoli, C., S. Zerbini, L. Lago, B. Richter, D. Simon,
F. Domenichini, C. Elmi, and M. Ghirotti. 2003. Influence of soil
consolidation and thermal expansion effects on height and
gravity variations. Journal of Geodynamics 35(4–5): 521–539.
DOI:10.1016/S0264-3707(03)00012-7.
[24] Seeber, G. 2003. Satellite geodesy, 2nd completely revised
and extended edition. Walter de Gruyter GmbH & Co. KG
10785: 303–304.
[25] Soler, T. and M. Chin 1985. On transformation of covariance
matrices between local cartesian coordinate systems and
commutative diagrams, In ASP-ACSM Convention, 393–406.
[26] U.S. Geological Survey. 2017. Earthquake lists, Maps, and
Statistics, accessed december 18, 2017 at.
[27] Wanninger, L. 2009. Correction of apparent position shifts
caused by GNSS antenna changes. GPS solutions 13(2):
133–139. DOI:10.1007/s10291-008-0106-z.
[28] Williams, S. 2003. The effect of coloured noise on
the uncertainties of rates estimated from geodetic
time series. Journal of Geodesy 76(9–10): 483–494.
DOI:10.1007/s00190-002-0283-4.
[29] Williams, S.D. 2008. CATS: GPS coordinate time series
analysis software. GPS solutions 12(2): 147–153.
DOI:10.1007/s10291-007-0086-4.
[30] Williams, S.D., Y. Bock, P. Fang, P. Jamason, R.M. Nikolaidis,
L. Prawirodirdjo, M. Miller, and D.J. Johnson. 2004.
Error analysis of continuous GPS position time series.
Journal of Geophysical Research: Solid Earth 109(B3).
DOI:10.1029/2003JB002741.
[31] World Data.info. 2021. The climate in iraq.
[32] Yenter, M., T. Wheeler, J. Mejia, K. Joyce, and M. Cline. 2005,
04. Developing IGRS. Engineer 36: 17–18.
[33] Zhang, J., Y. Bock, H. Johnson, P. Fang, S. Williams, J. Genrich,
S. Wdowinski, and J. Behr. 1997. Southern california
permanent GPS geodetic array: Error analysis of daily position
estimates and site velocities. Journal of geophysical research:
solid earth 102(B8): 18035–18055. DOI:10.1029/97JB01380.
