sensors
Article
Direction of Arrival Estimation of GPS Narrowband
Jammers Using High-Resolution Techniques
Mohamed Moussa 1,2 , Abdalla Osman 3 , Mohamed Tamazin 2, * , Michael J. Korenberg 1 and
Aboelmagd Noureldin 3
1
2
3
*
Electrical and Computer Engineering Department, Queen’s University, Kingston, ON K7L 3N6, Canada;
mohamed.moussa@queensu.ca (M.M.); korenber@queensu.ca (M.J.K.)
Electronics and Communications Engineering Department, Arab Academy for Science, Technology and
Maritime Transport (AASTMT), Alexandria, Egypt
Electrical and Computer Engineering, Royal Military College of Canada, Kingston, ON K7K 7B4, Canada;
osman.a.m7@gmail.com (A.O.); Aboelmagd.Noureldin@rmc.ca (A.N.)
Correspondence: tamazin@aast.edu
Received: 7 November 2019; Accepted: 11 December 2019; Published: 14 December 2019
Abstract: GPS jamming is a considerable threat to applications that rely on GPS position, velocity,
and time. Jamming detection is the first step in the mitigation process. The direction of arrival
(DOA) estimation of jamming signals is affected by resolution. In the presence of multiple jamming
sources whose spatial separation is very narrow, an incorrect number of jammers can be detected.
Consequently, mitigation will be affected. The ultimate objective of this research is to enhance
GPS receivers’ anti-jamming abilities. This research proposes an enhancement to the anti-jamming
detection ability of GPS receivers that are equipped with a uniform linear array (ULA) and uniform
circular array (UCA). The proposed array processing method utilizes fast orthogonal search (FOS)
to target the accurate detection of the DOA of both single and multiple in-band CW jammers. Its
performance is compared to the classical method and MUSIC. GPS signals obtained from a Spirent
GSS6700 simulator and CW jamming signals were used. The proposed method produces a threefold
advantage, higher accuracy DOA estimates, amplitudes, and a correct number of jammers. Therefore,
the anti-jamming process can be significantly improved by limiting the erroneous spatial attenuation
of GPS signals arriving from an angle close to the jammer.
Keywords: GPS; jamming detection; array processing; fast orthogonal search
1. Introduction
A GPS signal is at a disadvantage because of its low signal power, which is typically between
–125 dBm and –130 dBm at the earth’s surface and is therefore highly susceptible to intentional and
non-intentional interferences. GPS is relied upon in many critical infrastructure sectors [1]. Such an
example is the civil aviation sector, which is characterized by a safety-critical nature that imposes
stringent requirements on the availability of the GPS solution. In this case, GPS is used to help land
aircraft using ground-based augmentation systems (GBAS). They are steadily deployed, but their
operation can be disrupted when interference and jamming signals are in an airport’s vicinity. A clear
example of the threat is Newark Liberty International Airport, where the GBAS system’s operation was
affected due to GPS jamming using an illegal personal privacy device (PPD) [2]. When interference is
present, detecting it is a crucial first step toward its mitigation.
In order to detect and mitigate interference, many authors have proposed solutions that exploit
the time domain, frequency domain, space domains, or combinations of the fields. In the time
domain, pulse blanking is commonly used, whereas, in the frequency domain, adaptive filtering is
Sensors 2019, 19, 5532; doi:10.3390/s19245532
www.mdpi.com/journal/sensors
Sensors 2019, 19, 5532
2 of 18
used [3–5]. Time-domain pulse blanking and frequency excision introduce a set of drawbacks such
as distorted correlation peaks, acquisition problems, and tracking performance degradation. When
spatial mitigation is used, the drawbacks mentioned above will no longer be evident.
In the space domain, adaptive antenna arrays (smart antennas) that employ spatial processing
are used. Accordingly, several advantageous and unique capabilities are introduced, for example, the
direction of arrival (DOA) estimation, beam-steering, and null-steering. Once the DOA estimates are
computed, interference mitigation can be achieved using null steering by placing nulls in the detected
direction of interference [4,6–8]. On the other hand, beam-steering directs the main lobe toward the
direction of arrival of the GPS signal [9].
Therefore, smart antennas continue to attract interest in the field of GPS interference detection and
mitigation due to their superiority compared to single-antenna receivers, when dealing with intentional
and non-intentional interferences. The interference mitigation capability of smart antennas is mainly
based on the accuracy of the DOA estimation of interference sources [4,10,11]. For reliable operation,
GPS interference mitigation using smart antennas should be capable of detecting multiple interference
sources arriving at different DOAs. Conventional beam-forming methods provide DOA estimation for
interference sources with limited resolution. This degrades the spatial filtering capability of smart
antennas when the jamming signal is closely spaced to the DOA of the GPS signal. Accordingly,
high-resolution DOA estimation methods are more suitable for interference mitigation using smart
antennas. Various high-resolution DOA estimation methods have been developed, including MUltiple
SIgnal Classification (MUSIC) [12,13] and estimation of signal parameters via rotational invariance
techniques (ESPRIT) [14]. However, the performance of these methods degrades severely as the
coherence between jamming signals increases [15].
This paper aims to propose a new high-resolution DOA estimation method of jammers in order to
provide smart antennas with the most accurate estimate of the jammers’ DoAs such that beam-steering
and null-steering processors nullify the jammer signals without erroneously nullifying parts of the
GPS signal, thereby decreasing its effectiveness. The proposed method was designed to operate in
the precorrelation section of GPS receivers, where mitigation is more effective. It is based on the
orthogonal search (OS) technique [16] and creates bases for the steering function. The coefficients of the
steering function are determined via a unique orthogonal search method that provides high-accuracy
DOA estimation. The proposed method can enhance anti-jamming techniques that utilize antenna
arrays, specifically in situations where the DOA of jamming signals is very close to the DOA of the
GPS signal. Hence, it can enhance the null-steering process and avoid a reduction in GPS signal levels
while rejecting jamming signals. Furthermore, the high-resolution capabilities of the proposed method
can resolve jamming signals arriving at very close DOAs and hence improves the overall performance
of the anti-jamming method.
The proposed method is compared to classical DOA estimation and MUSIC using GPS L1 signals
obtained from a Spirent 6700 simulator; the arrays utilized are uniform linear array (ULA) and uniform
circular array (UCA) arrays. The jamming signals were simulated as CW sources originating from
different directions. The performance of the investigated methods is evaluated in the case of single and
multiple interference signals. Our results showed that the proposed method introduces a significant
enhancement to GPS jamming detection and hence can improve overall system performance.
2. DOA Estimation Methods
Spatial domain detection and mitigation is a mature field where many DOA estimation methods
have been proposed. Classical beam-forming, also known as conventional beam-forming, [17]
represents an earlier approach to DOA estimation; the spatial spectrum is scanned using predetermined
angular steps in search of the direction angle that corresponds to the highest spectral power. The
drawback of this method is its inability to form sharp peaks unless a large number of array elements
are used. Consequently, the method suffers from a limited ability to resolve closely spaced sources [18].
Sensors 2019, 19, 5532
3 of 18
Capon’s beamformer, also known as the minimum variance distortionless response (MVDR),
attempts to improve the drawbacks associated with the classical method and yield a significant
improvement [19]. The output power is minimized while constraining the gain in the necessary
direction to unity. The MVDR method requires an additional matrix inversion compared to the classical
method. Fortunately, it outperforms the conventional method in most cases. The disadvantages of the
MVDR method are the necessary other computations and the failure to estimate the DOAs of highly
correlated signals [20].
Advances in the field of direction-finding lead to the development of high-resolution DOA
estimation methods that utilize the signal subspace. By performing Eigenanalysis on the spatial
covariance matrix, the signal and noise subspaces are generated. MUSIC is among these subspace
techniques [21]. MUSIC showed that the steering vectors associated with the received signals are found
in the signal subspace. A search throughout the possible steering vectors is conducted; the ones that
are orthogonal to the noise subspace are designated as desired signals. However, errors arise in real
scenarios, since full orthogonality is challenging to achieve due to errors in the estimation of the noise
subspace. The MUSIC spectrum generates a tremendous value when a match between the generated
steering vector and the actual DOA occurs. The disadvantage of MUSIC is that it is not able to identify
DOAs of correlated signals on its own; the received signal must undergo a preprocessing technique
called spatial smoothing [22]. In addition, it is computationally expensive, sensitive to noise, and must
have prior knowledge of the number of signals it is looking for. Several variants of the method have
been proposed to improve its performance [23].
The estimation of signal parameters via rotational invariance technique (ESPRIT) algorithm
was proposed by [24]. It is computationally efficient and more robust if compared to MUSIC as it
does not search the entire spectrum, but a significant drawback is an incompatibility with all array
geometries, as it was designed for uniform linear arrays (ULA). The method has been extended to
include multidimensional arrays, and many versions have been produced [24,25].
Generally, spatial signal processing using antenna arrays is considered one of the most effective
techniques for narrowband interference detection and primarily suppression [7,8]. Antenna arrays
enable interference mitigation using null steering, beam steering, or a combination of both. This
is achieved using adaptive beam-forming and high-resolution DOA estimation methods [26]. The
obtained DOA estimates are utilized to produce nulls in the direction of interfering signals, and if
beam steering is available, steer the main beam toward the desired GPS signal.
Various DOA estimation algorithms have been developed for array signal processing applications.
The selection of the DOA estimation algorithm is a crucial element of adaptive antenna array design. It
directly affects null steering and beam-forming processes.
The DOA estimation techniques that are predominantly used are those that are based on classical
DOA, Capon, and MUSIC algorithms [26]. It has been reported that Eigen-decomposition based
methods such as MUSIC have high-resolution DOA estimation performance [26]. These methods
are based on exploiting the Eigenstructure of the input covariance matrix. Applying MUSIC to GPS
anti-jamming has shown to provide a significant enhancement to the overall system performance [27].
However, the performance of MUSIC is limited by the coherence of the jamming sources [28].
3. DOA Estimation Using Fast Orthogonal Search (FOS)
This research studies the application of FOS in the DOA estimation of jamming signals [5]. FOS is
a highly efficient general-purpose modeling method that has been used in several applications [29].
This method is a modification [15] of an original OS algorithm, which constructs a functional expansion
of an input data using an arbitrary set of non-orthogonal candidate functions. The functional expansion
of the input signal in terms of the arbitrary candidate functions is given by:
y(n) =
M
X
m=0
am pm (n) + e(n)
(1)
Sensors 2019, 19, 5532
4 of 18
where n is the sample index, {am }M
m = 0 are the weights of the functional expansion, and e(n) is the
modeling error. The signal arriving at the antenna array is generally composed of complex sinusoidal
signals and noise. The signal received at the nth sensor at time ti from S sources is defined by [30]:
ym (ti ) =
SX
−1
Pk (ti ) exp(j2πfk (ti − τm (θk ))) + e(ti )
(2)
k=0
where i = 0, . . . , T − 1 is the signal time index, m = 0, . . . , M − 1 is the sensor index number,
k = 0, · · · , S − 1 is the signal index number, θk is the DOA of the signal from kth source, Pk is the
complex amplitude of the kth source Pk (ti ) = Pk (ti ) exp(jϕ(ti )), ϕ(ti ) is the phase of the kth signal
at its source, and e (ti ) is the modeling error. For L real signals, S = 2 L and the complex amplitudes
are conjugate pairs (i.e., Pk + 1 = P∗k and fk + 1 = −fk ). The reference is taken at the first element;
thus, τ0 (θk ) = 0.
The following model can represent the received signal at the antenna array:
y(n) =
2L−1
X
ak (θk )sk (n) + e(n)
(3)
k=0
where ak (θk ) is the M × 1 antenna array steering vector for the kth source and sk is the kth source
signal received by the antenna array. For real signals, each candidate is represented by its amplitude,
sinusoidal function (sin or cos), and phase. In order to estimate the amplitude and phase of each
candidate, the fitting process is carried by selecting pairs of sin and cos functions for DOAs. Hence, for
L DOAs, we have 2L-1 candidate functions.
OS creates an orthogonal basis set [31] for the steering vectors using the Gram–Schmidt procedure.
This set forms a model equivalent to Equation (3), which is defined by:
y(n) =
2L−1
X
wk (θ0 , θ1 , · · · , θr )gk (n) + e(n)
(4)
k=0
where r = k/2 for k even values, r = (k − 1)/2 for k odd values, w∗i wj = 0 for i , j (* denotes
n
o
complex-conjugate transpose), and gk (n) are the coefficients of the orthogonal basis vectors.
The set of orthogonal basis vectors {wk } is obtained from the candidate functions ak (θk ) according
to the following formula:
k−1
X
wk = a0 −
αki wi
(5)
i=0
wi ∗ ak
wi ∗ wi
and w0 = a0 .
where αki =
The coefficients gk (n) are estimated using the following formula:
gk (n) =
wk ∗ y ( n )
.
wk ∗ wk
(6)
This formula was derived to minimize the mean square error (MSE) between the observations
vector y (n) and the orthogonal model. The total MSE (ε2 ) is defined as:
ε2 =
S−1
T−1
X
1X 1
gk (n)wk k2
k y(n) −
T
M
n=0
=
k=0
S−1
T−1
X
1 X ∗
Qk
y (n)y(n) −
MT
n=0
k=0
(7)
Sensors 2019, 19, 5532
5 of 18
where Qk represents the reduction introduced to the MSE by adding the candidate basis vector ŵk to
the model and defined by the following equation:
2
1 XT−1 [ wk ∗ y(n)]
wk H wk XT−1
Qk =
=
g (n).
n=0
n=0 k
MT
wk ∗ wk
MT
(8)
Equations (5) and (8) indicate that the calculation of Qk involves the correlations of wk with
themselves, the steering vectors ak , and the data y(n). The calculation of Qk involves the correlations
of wk with themselves; the steering vectors ak and y(n) are the dates. Therefore, the orthogonal search
can discard the need to create the orthogonal basis vectors wk explicitly.
This observation was used in [29,32,33] to accelerate the OS algorithm by discarding the explicit
creation of the orthogonal basis vectors wk . Accordingly, this approach was named FOS, which
stands for fast orthogonal search [29,32,33]. FOS calculates the correlations mentioned above using
Cholesky factorization. It starts by creating a variable Dkm , which represents the correlation between
the candidate and the steering vector. Thus, Dkm is defined by:
Dkm = wm ak = a∗m ak −
m−1
X
αmj ∗ Dkj
(9)
j=0
Dkj
where m = 1 . . . k and αkj = D .
jj
The correlation between the candidate basis vector wk and the observations vector y (n) is denoted
by Cnk and defined by:
k−1
X
(10)
αkj ∗ Cnj .
Cnk = wk ∗ y(n) = ak ∗ y(n) −
j=0
Consequently, the coefficients of orthogonal basis vectors gk (n) and the amount of reduction in
MSE Qk are expressed by:
C
gk (n) = nk
(11)
Dkk
Qk =
XT−1
n=0
Cnk 2
.
MTDkk
(12)
The derivation of Equations (11) and (12) are based on the properties of orthogonal basis vectors
wk , which imply that w∗k wk = w∗k ak . Additionally, given that w0 = a0 , the initial values for Dkm and
Cnk are defined by the following:
Dm0 = a∗0 am and Cn0 = a∗0 y(n).
(13)
FOS constructs the model defined by Equation (4) by adding the candidate functions corresponding
to DOA one at a time. The selection of candidate functions corresponding to an estimate θ̂k is based
on the assessment of a set of parameters to provide the value that yields maximum MSE reduction
Qk . Since a real signal is formed from accumulated sine and cosine components, the MSE reduction
contributed by a single candidate parameter can be represented by the sum of the MSE reductions
induced by the two corresponding steering vectors of that candidate (i.e., Qk and Qk+1 ). Thus, the
selected estimate θ̂k is defined by:
θ̂k = arg[max(Qk + Qk+1 )].
(14)
Once θk is selected, the two terms gk (n)wk and gk+1 (n)wk+1 are added to the constructed model.
This process is repeated until L parameter estimates are chosen. The values of the signal coefficients
sk (n) are obtained recursively using the following formula:
Sensors 2019, 19, 5532
6 of 18
ŝk (n) =
S−1
X
gm (n)vm
(15)
m=k
where vk = 1 and vm = −
m−1
P
j=k
αmj vj , m = k + 1, k + 2, . . . , S − 1.
3.1. Stopping Conditions for the FOS Algorithm
The FOS algorithm can be stopped using one of the following criteria:
a.
b.
c.
Reaching a predefined maximum number of terms to be fitted. This requires the knowledge of
the number of narrowband interference signals impinging upon the antenna array. However,
the orthogonal search approaches may be combined with any of several statistical criteria to
determine when to stop adding terms to the model, thus providing an estimate of the number of
signals [30].
Reaching a predefined threshold for MSE reduction. This criterion is achieved when the ratio of
MSE to the mean squared value of the input signal is below a predefined threshold. Accordingly,
the limitation of knowing the number of expected interference signals is waived. However, it
may lead to an increase in processing time. Thus, this ratio should be carefully selected to avoid
excessive processing time.
When adding another term to the model, the MSE reduction is less than the MSE reduction that
would be gained if white Gaussian noise was added.
3.2. FOS Complexity
FOS necessitates floating-point operations of the order of CMT + CL2 , where C is the number of
candidate steering vectors searched [33]. Moreover, if the elements of the array (and therefore, the
data samples) are not equally spaced, FOS will require a higher order, which becomes CMT + CL2 +
CML [30].
Candidate Function Selection for the FOS Algorithm
The process of candidate function selection plays a crucial role in the FOS algorithm. In this
research, the selected candidate functions are pairs of sine and cosine functions corresponding to the
DOAs of the search domain. Accordingly, the chosen candidate functions for a ULA and UCA using
the model shown in Equation (3) are given by [33]:
ak (θk ) = ak (m, θk )
(16)
where θk is the DOA of the kth source signal.
When θk is measured from the line of the antenna array, the following candidate functions are used
d
d
ak (m, θk ) = cos 2mπ
cos θk ak+1 (m, θk ) = sin 2mπ
cos θk .
λ
λ
(17)
On the other hand, when θk is measured from the line perpendicular to the antenna array, Equation
(18), which is shown below, is utilized:
d
d
ak (m, θk ) = cos 2mπ
sin(θk ak+1 (m, θk ) = sin 2mπ
sin(θk )
λ
λ
(18)
where m = 0, 1 . . . , M − 1, i = 0, 1 . . . , T − 1d is the spacing between elements, and λ is the wavelength
of the received signal.
Sensors 2019, 19, 5532
7 of 18
The coefficients sk (n) are given by [33]:
sk =
p
2Pr cos(πti + ϕr (ti ))sk+1 =
p
2Pr sin(πti + ϕr (ti ))
(19)
where Pr and ϕr (ti ) are the power and phase of the rth source, respectively.
4. Antenna Array Design
In the following sections, the array design for the ULA and UCA arrays will be studied from a
perspective that will ultimately serve the goal of DOA estimation. Therefore, in order to adhere to the
common literature notations, the total number of elements in an array will be referred to as N elements,
while n will represent the index of an element in the array and n = 0, 1, . . . , N − 1.
4.1. Uniform Linear Array
The ULA is capable of filtering the electromagnetic environment in which it operates, based on
the location of the signals in its vicinity. Array geometry has a significant effect on the DOA estimation
abilities of antenna arrays. In this study [5], the ULA shown in Figure 1 was considered.
M
M-1
dn
dn
dn
dn
dn
dn
Figure 1. Uniform linear array.
For the general ULA showed in Figure 1 with inter-element spacing and
dn dn and M identical
elements, the signal received by the antenna array is given by:
s(t)
.
M−1
X
−jk.d
n
e
y(t) = s(t)
y(t) = s(t)
e
(20)
n=0
where s(t) is the transmitted signal.
The quantity in parenthesis is known as the array factor (AF). This factoring is regularly termed
pattern multiplication. It can be utilized when the identical array elements are oriented in the same
direction. Moreover, the radiation pattern of the array is the result of the
of the element
dn multiplication
= (0,0, z ).
radiation pattern and array factor.
When the elements of the ULA are all located on the z-axis, dn = (0, 0, zn ). Then, the AF becomes:
AF =
n = 0, … , M − 1, and θ
e
M−1
X
2π
e−j λ zn cos θ
AF =
dn = (0,0,
AF = ∑
) n=0
(21)
e
where n = 0, . . . , M − 1, and θ is the direction
of arrival of the received signal.
For an M-element, ULA let dn = 0, 0, nλ
2 . Hence, the AF becomes:
∑
AF =
AF =
AF =
M
c =
XM−1
n=0
e−jnπ cos θ .
(22)
Sensors 2019, 19, 5532
8 of 18
The array factor formula can be simplified using the following identity [34]
1 − cM
.
1−c
(23)
1 − e−jMπ cos θ
.
1 − e−jπ cos θ
(24)
XM−1
n=0
cn =
Therefore, the array factor for ULA is given by:
AF =
After factoring, Equation (24) reduces to:
Mπ
Mπ
e−j 2 cos θ sin 2 cos θ
.
AF =
π
e−j 2 cos θ sin π2 cos θ
(25)
The AF magnitude pattern for ULA is shown in Figure 2. The magnitude of the array factor is
plotted for an array with M = nine elements. Based on the magnitude of the array factor, the maximum
energy is received or transmitted by the array when θ = 90◦ or when θ = 270◦ . Introducing weights and
varying the orientation of the ULA can manipulate this array factor. Since the array factor is a linear
function of the weights, choosing weights is considered a minor procedure. On the other hand, the
array factor is a nonlinear exponential function of the element positions; consequently, array geometry
optimization is more challenging [34].
90
10
120
60
8
6
150
30
4
2
180
0
210
330
240
300
270
Figure 2. Array factor for uniform linear array (ULA).
Other vital parameters of array factors include the beam width and sidelobe level. The beam
width is usually defined as a null-to-null or half-power beam width. The null-to-null beam width is
the angular distance between the first nulls around the main beam [34]. The half-power beam width is
the angular distance between the half-power points (i.e., 3 dB points on the array factor) around the
main beam. The sidelobe level is usually defined as the maximum value of the array factor that is
found outside the main beam.
4.2. Uniform Circular Array
The planar array geometries can be divided into three other subcategories: circular, rectangular,
and square. The circular arrays stand out of the three, as they do not have edge elements. In the
absence of edge constraints, the beam pattern of a circular array can be electronically rotated. Moreover,
the circular arrays also have the capability to compensate for the effect of mutual coupling by breaking
down the array excitation into a series of symmetrical spatial components [35,36].
The array factor of a circular antenna array, which is shown in Figure 3, is given by [35]:
𝜙
AF(𝜃, 𝜙) =
𝑘
e(
(
𝑛 𝑀
𝑘=
))
𝑏
Sensors 2019, 19, 5532
9 of 18
AF(θ, φ) =
M
X
ej(kb sin θ cos (φ−φn ))
(26)
n=1
where φn is the angular position of element n, M is the number of antenna elements, b is the radius of
the circle, and k is the wavenumber, k = 2π
λ .
Figure 3. Circular array geometry.
The array factor for a circular antenna array is a function of both spherical angles; therefore, it can
filter the received signals based on their azimuth and elevation angles. The effect of the array on the
received signal as a function of the angle of arrival can be demonstrated by examining the array factor.
The magnitude pattern of the circular array factor is shown in Figures 4–6. This pattern is drawn for a
circular array of nine elements that are equally spaced on a circular circumference with a radius of
λ
2 . The circular arrays configuration does not have the edge effect. It is clear from the AF magnitude
shown in Figure 6 that the circular array response is symmetrical around the elevation direction, which
enhances the capabilities of circular antenna arrays in beam steering.
Figure 4. Array factor magnitude for the circular array.
Sensors 2019, 19, 5532
10 of 18
Figure 5. Array factor pattern for the circular array.
Figure 6. XY-plane view of an array factor pattern for the circular array.
5. Results and Discussion
The following section demonstrates the results obtained from the different scenarios that were
implemented in order to evaluate the performance of the developed high-resolution DOA estimation
method for ULA and UCA arrays geometry. The proposed method was tested with single and multiple
jamming signals at different jamming-to-signal ratios (JSRs). The GPS L1 signals were obtained from a
Spirent GSS6700 simulator and were down-converted and digitized using the Novatel Firehose digital
frontend, as shown in Figure 7.
Sensors 2019, 19, 5532
11 of 18
Spirent GSS6700 Multi-GNSS
Simulator
LNA
SimGENTM
Software
NovAtel
FireHose
Front-end
MATLAB
Figure 7. Experimental setup of the Spirent GSS6700 and the Firehose.
LNA
Referring to [37], one radio frequency (RF) front end and a single RF output Spirent GSS6700
simulator running SimGENTM software are utilized to verify the proposed method. Initially, the
desired experimental parameters such as the number of antenna elements, array geometry, GPS
signal frequency (i.e., L1, L2 and L5), and consequently antenna element spacing are defined. Upon
designating the zero-phase element’s (reference element) location in terms of latitude, longitude, and
height, the remaining elements of the array are mapped according to the predefined experimental
configuration. Once
the locations of all the antenna elements are set, a simulation scenario is created
SimGENTM
using theSoftware
simulator software. The simulator software enables the reiteration of the simulation scenario
without changing any parameters except for the location of the antenna element in 3D space. Figure 8
illustrates the simulation sequence for a given M-element array, where each iteration corresponds to
one of the array elements where n = 1, 2 . . . , M.
Figure 8. Simulation sequence of an M-element array of GPS antennas [37].
The jamming signals frequencies are relative to the GPS signal’s center frequency obtained from
digitization and down-conversion. Jamming signals were simulated by adding sinusoidal signals to the
DGA output. The simulated jammer frequencies ranged from 100 to 400 Hz. This range ensures that the
introduced jamming signals fall within the bandwidth of the DGA output digitized, down-converted
GPS L1 signals, which are baseband signals with 1.023 MHz bandwidth.
5.1. Results for ULA
•
Single Jammer
𝑀
𝑛 = 1,2 … , 𝑀
The results obtained for a single jammer simulated at a frequency of 100 Hz and arriving at an elevation
angle of 50◦ at a JSR of 15 and 45 dB are shown in Figures 9 and 10, respectively. The figures illustrate
the normalized output for the classical, the proposed, and MUSIC DOA estimation methods. The
Sensors 2019, 19, 5532
12 of 18
figures demonstrate that all three methods detected the single jammer with high accuracy at JSRs of 15
and 45 dB. When the detection of a single GPS jamming signal is required, the performance of all three
methods is satisfactory, as the jammer of interest’s received power is usually high.
1
Classical
Proposed
Classical Method
MUSIC
Proposed Method
Normalized
Normalized
Spatial
Spatial
Spectrum
Spectrum
1
0.9
0.9
0.8
MUSIC
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
0
20
40
0
20
40
60
80
100
120
60 Angle
80 in Degrees
100
120
Angle in Degrees
140
160
180
140
160
180
Figure 9. Direction of arrival (DOA) estimation of one jammer at a JSR = 15 dB.
1
Classical
Proposed Method
Classical
MUSIC
Proposed Method
MUSIC
Normalized
Normalized
Spatial
Spatial
Spectrum
Spectrum
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
0
20
40
0
20
40
60
80
100
120
60 Angle
80 in Degrees
100
120
Angle in Degrees
140
160
180
140
160
180
Figure 10. DOA estimation of one jammer at JSR = 45 dB.
The main contribution of the proposed method when estimating the DOA of a single jammer is
the high accuracy in the detection of jamming signals’ amplitude, as shown in Figure 10. This is mainly
related to the nature of the operation of FOS on which the method is built. It operates by constructing
a signal model determined by candidate functions corresponding to the detected signals.
•
Multiple Jamming Signals
The performance analysis of FOS versus MUSIC and classical DOA was examined in terms of resolution,
tolerance to JSR, and tolerance to the jamming signals coherence. For this purpose, three jamming
signals were simulated with frequencies of 100 Hz, 200 Hz, and 300 Hz and arriving at elevations of
Sensors 2019, 19, 5532
◦
◦
13 of 18
50° , 60° , and 70°
◦
50 , 60 , and 70 , respectively. The tolerance to jamming signals coherency was examined by repeating
this scenario with three jamming signals at frequencies of 100, 105, and 11050
Hz.
°
, 60° , and 70°
◦
◦
The results obtained for three jammers arriving at elevations of (50 , 60 and 70◦ ) JSRs are shown
in Figures 11 and 12. The figures indicate the normalized output for classical DOS estimation, MUSIC,
and the proposed DOA estimation methods. The JSRs of the jamming signals used in Figures 11 and 12
are 15 and 45 dB, respectively.
1
Classical
Proposed Method
MUSIC
Normalized Spatial Spectrum
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
20
40
60
80
100
120
Angle in Degrees
140
160
180
Figure 11. DOA estimation of three jammers at JSR = 15 dB and frequencies of 100 Hz, 200 Hz, and
300 Hz.
1
Classical
Proposed Method
MUSIC
Normalized Spatial Spectrum
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
20
40
60
80
100
120
Angle in Degrees
140
160
180
Figure 12. DOA estimation of three jammers arriving at JSR = 45 dB and frequencies of 100 Hz, 200 Hz,
and 300 Hz.
The results shown in the above figures demonstrate that the proposed FOS-based method
outperformed both the MUSIC and classical methods due to its high tolerance to the variation of JSR
and jamming signals coherency. MUSIC detected three jammers at a JSR of 45 dB with zero error in
Sensors 2019, 19, 5532
14 of 18
estimated DOAs, but its performance degraded at a JSR of 15 dB as it detected two jammers only. On
the other hand, our method’s performance was stable in terms of the number of jammers detected at
JSRs of 15 dB and 45 dB, as it detected three jammers accurately with zero error in estimated DOAs.
It also degraded slightly at a JSR of 15 dB where its power allocation for detected jammers was not
equally divided among the three jammers that were simulated with equal power.
5.2. Results for UCA
The proposed method (FOS) performance was examined using a UCA configuration described
earlier, with seven elements equally spaced on the circumference of a circle with a radius of λ2 .
•
Single Jammer
A single jammer was simulated as a 100 Hz sinusoidal signal arriving at an elevation of 40◦ and an
azimuth of 40◦ . Figures 13–15 demonstrate that the performances of the FOS, MUSIC, and classical
DOA in 2D single jammer detection are almost identical and show a clear detection of the jamming
signal with zero error in both elevation and azimuth. The advantage of FOS and MUSIC is that their
spatial spectrum has much higher resolution compared to that of classical DoA.
Figure 13. DOA estimation of one jammer using classical DoA at JSR = 45 dB.
Figure 14. DOA estimation of a single jammer using MUSIC at JSR = 45 dB.
Sensors 2019, 19, 5532
15 of 18
Figure 15. DOA estimation of one jammer using fast orthogonal search (FOS) at JSR = 45 dB.
•
Multiple Jamming Signals
In this scenario, three equal power jammers were simulated at frequencies of 100, 200, and 300 Hz
with elevations and azimuths separations of 10◦ . The DoAs of the three jammers were simulated at
(30◦ , 30◦ ), (40◦ , 40◦ ), and (50◦ , 50◦ ). This experiment is conducted at JSRs of 15 and 45 dB. The results
obtained from the MUSIC and classical DoA estimation methods are shown in Figures 16 and 17. In
this scenario, the performance of FOS is compared to MUSIC only, as their performance versus the
classical DoA method was already examined for ULA and showed superior performance in all cases
over classical DoA.
Figure 16. DOA estimation of three closely spaced jammers using MUSIC at JSR = 15 dB.
The results shown in Figure 17 showed steady performance for FOS in terms of the resolution
and number of detected jammers for the cases where the JSR was 15 dB. The degradation in FOS was
observed in its allocation of power to jamming signals, as it was not accurately allocated for jammers
arriving at a JSR of 15 dB. The increase in the level of noise affected FOS performance and caused most
of the power to be allocated to the jammer arriving at 40◦ . MUSIC had much lower accuracy in power
allocation for the three jamming signals. The degradation in MUSIC performance was more evident at
lower JSR ratios, where MUSIC was able to detect only two jamming sources, as shown in Figure 16.
The accuracy of the detected DOA elevation and azimuth of jamming signals at 15 dB JSR are close for
both FOS and MUSIC, as shown in Table 1, with a slight increase of 1◦ in the error of jamming signals
detected by MUSIC.
Sensors 2019, 19, 5532
16 of 18
Figure 17. DOA estimation of three closely spaced jammers using FOS at JSR = 15 dB.
Table 1. Results of three closely spaced jammers at JSR = 15dB.
Jammer (s)
1
2
3
Elevation
Azimuth
MUSIC
MUSIC (Error)
FOS
FOS (Error)
30◦
40◦
30◦
40◦
(33◦ , 33◦ )
(3◦ , 3◦ )
(29◦ , 29◦ )
(1◦ , 1◦ )
N/A
N/A
(40◦ , 40◦ )
(0◦ , 0◦ )
50◦
50◦
(48◦ , 48◦ )
(2◦ , 2◦ )
(50◦ , 50◦ )
(0◦ , 0◦ )
6. Conclusions and Future Work
High-resolution DOA estimations for interference detection using the ULA and UCA
configurations were studied. A new procedure was utilized to simulate GPS received signals
using a single RF output Spirent GSS 6700 GPS simulator and a single RF input Novatel FireHose
frontend. Postprocessing validated the GPS signals used. The main contribution of this research is the
application of nonlinear signal modeling techniques for GPS jamming detection using ULAs and UCAs.
The method employed for jammer DOA estimation is based on FOS. Its performance was evaluated
in terms of DOA estimation accuracy, magnitude estimation, and the number of jammers detected
at different JSRs. The proposed method was compared to MUSIC and classical DOA estimation. It
demonstrated significant improvement in estimating the magnitude of the spatial spectrum when
compared to MUSIC and consistently outperformed it. MUSIC only showed similar performance when
more powerful jamming signals were applied. The proposed method was more accurate at detecting
the amplitudes of the single jammers. In the case of ULA, as for multiple jammers, it outperformed
MUSIC by 33% in terms of the number of jammers detected at relatively low JSR. In the case of UCA,
as for multiple jammers, the proposed method detected all three jammers, and MUSIC only detected
two at relatively low JSRs. Moreover, the proposed method also demonstrated higher accuracy in the
detection of the three sources of jammers with slight variations in their amplitude at a relatively low
JSR of 15 dB. In addition, the capabilities of the classical method are limited by its low resolution, and
MUSIC is much more sensitive to noise and performs worse at lower signal power. Since the proposed
method is not optimized for real-time operation, it is limited by its higher computational complexity
compared to the classical approach. Nevertheless, the added complexity is justified by the high gains
in accuracy. The proposed method will detect a higher number of closely spaced jammers and achieve
higher accuracies as the number of array elements increases.
Author Contributions: M.M. and A.O. proposed the idea. M.M. and M.T. conceived and designed the simulations,
analyzed the data, and wrote the paper. M.J.K. and A.N. reviewed the paper and simulation results.
Sensors 2019, 19, 5532
17 of 18
Funding: This research received no external funding.
Conflicts of Interest: The authors declare no conflict of interest.
References
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
Goward, D.A. Expert Advice: The Low Cost of Protecting America. GPS World, 1 January 2014.
Guenter, D.; Dennis, J. Initial operational experience with CAT I Ground Based Augmentation System
(GBAS). In Proceedings of the 2015 Integrated Communication, Navigation and Surveillance Conference
(ICNS), Herndon, VA, USA, 19–24 April 2015.
Ojeda, O.A.Y.; Grajal, J.; Lopez-Risueno, G. Analytical Performance of GNSS Receivers using Interference
Mitigation Techniques. IEEE Trans. Aerosp. Electron. Syst. 2013, 49, 885–906. [CrossRef]
Li, M.; Dempster, A.G.; Balaei, A.T.; Rizos, C.; Wang, F. Switchable Beam Steering/Null Steering Algorithm
for CW Interference Mitigation in GPS C/A Code Receivers. IEEE Trans. Aerosp. Electron. Syst. 2011, 47,
1564–1579. [CrossRef]
Moussa, M. High Resolution Jamming Detection in Global Navigation Satellite Systems. Ph.D. Thesis,
Queen’s University, Kingston, ON, Canada, 2015.
Li, Q.; Wang, W.; Xu, D.J.; Wang, X.P. A Robust Anti-Jamming Navigation Receiver with Antenna Array and
GPS/SINS. IEEE Commun. Lett. 2014, 18, 467–470. [CrossRef]
Graas, F.V.; Soloviev, A. Beam Steering in Global Positioning System Receivers using Synthetic Phased
Arrays. IEEE Trans. Aerosp. Electron. Syst. 2010, 46, 1513–1522.
Javier, A.; Carles, F.-P.; Pau, C. Antenna Array Based GNSS Signal Acquisition for Interference Mitigation.
IEEE Trans. Aerosp. Electron. Syst. 2013, 49, 223–243.
Kaplan, E.D.; Hegarty, C.J. Understanding GPS/GNSS: Principles and Applications; Artech House: Norwood,
MA, USA, 2017. [CrossRef]
Trees, H.L. Optimum Array Processing Part IV of Detection Estimation, and Modulation Theory; John Wiley and
Sons, Inc.: New York, NY, USA, 2002.
Dempster, A.G.; Cetin, E. Interference Localization for Satellite Navigation Systems. Proc. IEEE 2016, 104,
1318–1326. [CrossRef]
Maloney, J.A.; Kwon, D.H.; Keller, S.D.; Janaswamy, R. Realistic GPS Coverage Prediction for Dual-Polarized
Controlled Reception Pattern Antennas. IEEE Antennas Wirel. Propag. Lett. 2017, 16, 19071910. [CrossRef]
Lin, T.; Abdizadeh, M.; Broumandan, A.; Wang, D.; O’Keefe, K.; Lachapelle, G. Interference Suppression
for High Precision Navigation Using Vector-Based GNSS Software Receivers. In Proceedings of the 24th
International Technical Meeting of the Satellite Division of the Institute of Navigation, Portland, OR, USA,
20–23 September 2011; pp. 372–383.
Mewes, H. Practical Comparisons Between Various High Resolution Spectral Estimation Algorithms for
Radio Direction Finding. In Proceedings of the 1995 Ninth International Conference on Antennas and
Propagation, ICAP ’95, Eindhoven, The Netherlands, 4–7 April 1995.
Rieken, D.W.; Fuhrmann, D.R. Generalizing MUSIC and MVDR for Multiple Noncoherent Arrays. IEEE
Trans. Signal Process. 2004, 52, 2396–2406. [CrossRef]
Korenberg, M.J. A robust orthogonal algorithm for system identification and time-series analysis. Biol.
Cybern. 1989, 60, 267–276. [CrossRef]
Rao, B.R. GPS/GNSS Antennas; Artech House: Norwood, MA, USA, 2013.
Liberti, J.C.; Rappaport, T.S. Smart Antennas for Wireless Communications: IS-95 and Third Generation CDMA
Applications; Prentice Hall PTR: Upper Saddle River, NJ, USA, 1999; p. 528.
Bellofiore, S. Smart Antenna Systems for Mobile Platforms. Ph.D. Thesis, Arizona State University, Tempe,
AZ, USA, 2002.
Johnson, D.H. The application of spectral estimation methods to bearing estimation problems. Proc. IEEE
1982, 70, 1018–1028. [CrossRef]
Schmidt, R. Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propag. 1986,
34, 276–280. [CrossRef]
Xu, W.; Kaveh, M. Comparative study of the biases of MUSIC-Like estimators. Signal Process. 1996, 50, 39–55.
[CrossRef]
Sensors 2019, 19, 5532
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
18 of 18
Barabell, A. Improving the resolution performance of eigenstructure-Based direction-Finding algorithms. In
Proceedings of the ICASSP ’83. IEEE International Conference on Acoustics, Speech, and Signal Processing,
Boston, MA, USA, 14–16 April 1983; pp. 336–339.
Roy, R.; Kailath, T. ESPRIT-estimation of signal parameters via rotational invariance techniques. IEEE Trans.
Acoust. Speech Signal Process 1989, 37, 984–995. [CrossRef]
Guanghan, X.; Silverstein, S.D.; Roy, R.H.; Kailath, T. Beamspace ESPRIT. IEEE Trans. Signal Process. 1994, 42,
349–356. [CrossRef]
Wang, Q.; Wang, L.; An, K.; Shou, Z.; Zhang, H. DOA Estimation of Smart Antenna Signal Based on MUSIC
Algorithm. J. Netw. 2014, 9, 1309–1316. [CrossRef]
Applebaum, S. Adaptive arrays. IEEE Trans. Antennas Propag. 1976, 24, 585–598. [CrossRef]
Vasyiyshyn, V.I. Antenna array signal processing with high-resolution by modified beamspace music
algorithm. In Proceedings of the 2007 6th International Conference on Antenna Theory and Techniques,
Sevastopol, Ukraine, 17–21 September 2007; pp. 455–457.
Korenberg, M.; Brenan, C.H.; Hunter, I. Raman spectral estimation via fast orthogonal search. Analyst 1997,
122, 879–882. [CrossRef]
Adeney, K.; Korenberg, M. Fast orthogonal search for array processing and spectrum estimation. IEE Proc.
Vis. Image Signal Process 1994, 141, 13–18. [CrossRef]
Humphreys, T.E.; Ledvina, B.M.; Psiaki, M.L.; O’Hanlon, B.W.; Kintner, P.M., Jr. Assessing the spoofing
threat: Development of a portable GPS civilian spoofer. In Proceedings of the ION GNSS international
technical meeting of the satellite division, Savannah, GA, USA, 16–19 September 2008; p. 56.
Korenberg, M.J. Identifying nonlinear difference equation and functional expansion representations: The fast
orthogonal algorithm. Ann. Biomed. Eng. 1988, 16, 123–142. [CrossRef]
Adeney, K.M.; Korenberg, M.J. Fast orthogonal search for direction finding. Electron. Lett. 1992, 28, 2268–2269.
[CrossRef]
Bevelacqua, P.J. Antenna Arrays: Performance Limits and Geometery Optimization; Arizona State University:
Tempe, AZ, USA, 2008.
Noordin, N.H.; Zuniga, V.; El-Rayis, A.O.; Haridas, N.; Erdogan, A.T.; Arslan, T. Uniform Circular Arrays for
Phased Array Antenna. In Proceedings of the 2011 Loughborough Antennas & Propagation Conference,
Loughborough, UK, 14–15 November 2011.
Dalli, A.; Zenkouar, L.; Adiba, E.F.; Bri, S. Circular Array with Central Element for Smart Antenna. Electr.
Electron. Eng. 2013, 3, 86–95.
Noureldin, A.; Moussa, M.M.E.; Osman, A.; Tamazin, M. Analysis of Different GNSS Array Processing
Methods Utilizing New Experimental Approach Using a Spirent Simulator and Single Frontend Receiver. In
Proceedings of the 2016 11th International Conference on Computer Engineering & Systems (ICCES), Cairo,
Egypt, 20–21 December 2016.
© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access
article distributed under the terms and conditions of the Creative Commons Attribution
(CC BY) license (http://creativecommons.org/licenses/by/4.0/).