An optimal design of GNSS interference localisation wireless security network based on time-difference of arrivals for the Arlanda international airport

1 Introduction

Global navigation satellite systems (GNSSs) are common tools for navigating different types of vehicles, including aircrafts. Receiving healthy signals, free of any interference, is a necessity for a successful navigation process. Navigation of aircrafts needs more attention as interference in their navigation signal may lead to catastrophes for crew, passengers and even people on the ground. Landing is a risky part of a flight around an airport, especially when the airport is not fully visible from the aircraft. Weather conditions, connection with the airport traffic control tower as well as healthy navigation signals are important factors for a successful landing. Any intentional or unintentional signal interference might lead to serious problems and risk people’s life. Establishment of wireless security networks of sensors, or GNSS receivers with the possibility of providing information about interference, over airports is a necessity today. From the information received by these sensors and their positions, the interference device can be localised with some uncertainties. The main issue is to select an optimal location for these sensors, which should be to achieve the best possible coverage. Here, a constrained quadratic optimisation model based on time-difference of arrivals (TDOA) is developed and applied for optimal estimation of the location of these sensors, or anchor nodes (ANs), in such a way that any interference device can be localised with higher precision or less level of uncertainties.

Two well-known types of signal interference are jamming and spoofing; the former deals with transmitting a signal into the same band as, or a band nearby to, the satellite navigation band of interest to jam it and the latter the transmission of a fake GNSS signal (Dempster 2016). There are studies showing that a simple and relatively cheap GNSS spoofer can be used to overtake, for example, a ship navigation without being detected (Humphreys et al. 2008 and Divis 2013). Since the power level of the GNSS signals is low, such signals are susceptible to interference; therefore, a relatively weak interference signal can jam a receiver (Dempster 2016). There are real examples that this interference affected operational infrastructures (Balaei et al. 2007, Clynch et al. 2003, Grant et al. 2009, Hambling 2011, Motella et al. 2008 and Pullen et al. 2012). Specifically, we can point to the unintentional cases range from a faulty TV amplifier, which jammed the Global Positioning System (GPS) operation at a harbour in Monterey, California, for 37 days (Clynch et al. 2003). A small jammer, which was used in a delivery van, disrupted the ground-based augmentation system aiding aircraft approaches at Newark Airport while driving on a nearby highway in 2009 (Hambling 2011, Pullen et al. 2012 and Warburton and Tedeschi 2011). The Central Radio Management Office of South Korea reported several disruptions from 2010–2012 due to GPS jammers being affected (Seo and Kim 2013). In Australia, Balaei et al. (2007) detected some interference and in Italy some from TV signals in the GNSS band, disrupting GPS (Motella et al. 2008). Recognition of an interference signal amongst all scattered signals is a complicated process and required skills in signal processing, which is outside the scope of this article.

Localisation of an interference device can be done from ANs equipped by GNSS receivers, which can detect interference and measure the signal time of arrival (TOA), angle of arrival (AOA) and TDOA from the device. By the known coordinates of the ANs and these measurements, the coordinates of the device are determined, or in other words, the device is localised. Drake and Dogancay (2004) performed this process by prolate spheroidal coordinates and stated that the mathematical equations of TDOA will be greatly simplified in the case of using these coordinates. However, as will be shown in this study, the range difference equation does not have complicated mathematical formula for estimating the interference device coordinates. Ananthasubramanian and Madlhow (2008) investigated AOA and developed a sequential algorithm and concluded that the localisation error is proportional to the AOA error variance, coverage area and reducible by increasing the number of estimates. According to the least-squares principle, when the redundancy of the system of equations increases the variance decreases, in addition, the error of localisation is always proportional to the error of localisation. Also, localisation using AOA needs antenna arrays (Trinkle et al. 2012) for mathematical modelling and estimation of AOA from these arrays (Huang et al. 2022). Thompson et al. (2009) studied the optimal configuration of the sensors' location and presented a method using differences-of-received-signal-strength measurements. They concluded that these measurements could be alternatives to TAO, AOA, and TDOA. Thompson (2013) investigated interference detection and localisation by analysing the dilution of precision (DOP), from received signal strength and TDOA and concluded that TDOA are superior to the received signal strength measurements.

Eshagh (2022) pointed out that localisation with TOA requires precise time synchronisation of transmitter and receivers so that distances can be computed from the measured TOAs. However, in a 2D localisation by using at least three sensors/receivers, the transmission time can be considered as an extra unknown in the system of equations and approximated simultaneously with the coordinates of the transmitter. The advantage of using TDOA is that the transmitter needs no synchronisation with the receivers/sensors (Gustafsson 2010, p. 78). The TDOAs between the ANs are estimated by cross-correlation processes amongst the received signals (Lindström et al. 2007). In addition, there are new GNSS environmental monitoring systems consisting of several low-cost sensors to monitor GNSS system performance in a specific area (Trinkle et al. 2012).

In the classical geodetic networks, an optimal configuration is obtained by maximising the precision and reliability of the network (Xu 1989, Koch 1982, 1985, Kuang 1996, Eshagh and Kiamehr 2007, or Eshagh and Alizadeh-Khameneh 2015). A wireless sensor network can be regarded also a type of geodetic network, but with different observables or structures. The control nodes (CNs) or points vary in such a way that the desired configuration is achieved in geodetic networks, whilst a wireless localisation security network, the ANs are displaced to reach to the optimal configuration. Eshagh (2022) developed a quadratic optimisation method based on TOA, AOA and TDOA with three-ANs and applied it at the Landvetter international airport in Sweden. However, his optimisation with TDOA over the airport was not successful because lack of enough number of nodes and the airport special geometric shape. TOA, AOA, and TDOA of signals to the sensors were considered observables, and a criterion matrix was selected for the precision of the CNs over the airport. The ANs’ coordinates varied until the estimated variance–covariance (VC) matrices of the CNs are fitted, but subjected to some constraints, to this pre-defined criterion matrix

This article is a continuation of the studies of Eshagh (2022), for the optimal design of configuration created by ANs, with the following differences:

• Four ANs are used for optimisation, which means to have a higher redundancy.

• More details about optimisation using TDOA are presented.

• Developing and applying directional constraints, which keeps the ANs along the airport runways, which has not been done in any optimisation study so far.

• Applying the developed method for the Arlanda international airport of Sweden.

2 GNSS interference localisation security network

A GNSS interference localisation wireless security network consists of a series of point covering a control area. These points are named CNs, and they are probable locations of the interference transmitter. They have known coordinates according to the resolution of the grid of these points or nodes in a pre-defined local coordinate system. ANs refer to the nodes with known coordinates and are supposed to have sensors/receivers. Figure 1 is a schematic GNSS interference localisation security network with small black circles as CNs and four ANs shown by triangles. The geometric form, the quadrilateral, created by the ANs is named configuration.

Figure 1

Interference localisation security network with four ANs.

3 TDOA and localisation from ANs

An interference signal reaches to ANs at different times based on the range difference from the ANs. By measuring these time-differences and the known signal speed, the range differences are obtained. The mathematical formula of a range difference (d _ijk ) between an interference device at the jth CN and ith and kth the ANs is

d i j k = L i j − L i k = ( x i − x j ) 2 + ( y i − y j ) 2 − ( x k − x j ) 2 + ( y k − y j ) 2 , j , k = M , N , O , P , (1a)where x _j and y _j are coordinates of the interfering device, and x _i , y _i and x _k , y _k are the pair coordinates of ith and kth ANs. L _ij and L _ik are, respectively, the distances of the ith and kth ANs to the interference device at the jth CN.

Since there are four ANs in our design, six mathematical formulae of TDOA can be constructed. The purpose is to estimate the coordinates of the interference device. However, due to the nonlinearity of the mathematical models, they ought to be linearised by the Taylor series around some approximate coordinates of the device. By linearising all these six models, the following system of equations of Gauss–Markov type is constructed for estimating the coordinates for the device (Koch 2010):

where ε is the vector of random errors with E{ε} = 0, where E{} stands for the statistical expectation, C _L the VC matrix of TDOA, Q the co-factor matrix carrying the geometrical properties of the network, and finally σ 0 2 a priori variance of unit weight. A is the coefficient matrix having partial derivatives of each formula of TDOA with respect to the coordinates of the device, for instance, in this study, the assumption is to measure six TDOA for estimating 2D coordinates of the interference device having two components of x and y, i.e. A is 6 × 2

where j = M, N, O, and P.

x is a vector of the coordinate updates to their initial values of x and y, L vector of differences between the observed and computed TDOA from the initial coordinates

where d i j k 0 , j, k = M, N, O, P is the range differences, and L i j 0 , j = M, N, O, P are the ranges computed based on the approximate coordinates of the ith CN.

The least-squares solution of equation (1b) is (Cooper 1987):

and the VC matrix of the estimated coordinates:

with the a posteriori variance of unit weight:

Since six TDOA can be used and the unknowns are the 2D coordinates of an interference device, the redundancy of the system will be four.

In the case where the TDOAs are measured from the interference device and to each AN pair, the coordinates of the device and their errors can be simply estimated. However, the main objective of this study is optimisation of the configuration of the ANs in such a way that the device can be located with higher precision and not localisation. This means that for the design, no measurement of TDOA and no approximate value for the interference device are needed and the matrix A, carrying the geometrical properties of the network, would be enough as the coordinates of the probable location of it are already determined from the grid of CNs.

4 Optimal configuration of ANs

Equation (1f) is the VC matrix of coordinates of the interference device, which is dependent on the coordinates of the ANs. Now, assume that the coordinates of the CNs and ANs are known from the grid. Therefore, computing the initial VC matrices of all CNs will be straightforward. After that, the ANs’ configuration is optimised by varying their coordinates during the optimisation process.

The initial VC matrix of each CN is 2 × 2. A diagonal matrix with the same dimension is normally used as a criterion for the precision of this point, that is equal variances for x and y coordinates and no covariance between them. All initial VC matrices of CNs should be fitted to this criterion by varying the coordinates of the ANs in a least-squares sense. Let us first linearise this VC matrix by the Taylor series:

where C x ˆ is the criterion matrix and C x 0 is the initial VC matrix derived from the initial coordinates of the ANs and a CN, Δx _j and Δy _j are the coordinate updates for optimising ANs’ configuration. The partial derivatives of C x 0 with respect to the x- and y-coordinates of the jth AN have the following expression:

where

The structures of derivatives of A are:

where

From equation (2a), one can conclude that Δx _j , Δy _j should be estimated in such a way that C x 0 , computed from the update coordinates, is fitted to the desired C x ˆ . To do so, for each element of the VC matrix, an equation should be constructed, and since total number of these elements is four for each CN, then the constructed system of equations will have four rows and since the number of ANs is four and each one has two coordinates, then the number of unknown parameters will be eight. Consequently, our system of equations is underdetermined. The number of ANs is constant in the whole wireless security network; therefore, by adding the elements of the VC matrix of another CN, the system will have four new additional equations, and finally by using all CNs, the number of equations will be four times of the number of CNs and much higher than eight and the system becomes overdetermined.

This overdetermined system is presented in the following form:

where ε′ vector of residuals, and ΔL is the vector of differences between the elements of the criterion and initial VC matrices, B stands for the coefficients matrix containing the partial derivatives of the VC matrix with respect to the ANs’ coordinates, and Δx stands for the ANs’ coordinate updates. The mathematical descriptions of them are:

where operator “vec” insert the columns of the VC matrices below each other and convert the 2 × 2 matrices to 4 × 1 vectors and ()^T stands for transposition operator of matrix algebra, n means the number of the control points.

The objective function for such an optimisation problem will be:

This minimisation of equation (3e) leads to the least-squares solution for Δx. However, the issue is that there is no control over the estimation of the coordinates of ANs, and they may move towards each other and make the system of equation (3a) ill-conditioned, or they may move outside that study area or become colinear. Consequently, this minimisation problem requires some constraints for controlling the movements of the ANs during the iterative optimisation process.

5 Limiting search area of ANs

Liming search areas around the initial positions of ANs is an effective way to control the movements of the ANs. A search area around the jth AN is defined by

where w j L and w j U are, respectively, the lower and upper bounds of the inequality constraints; and v j L and v j U are corresponding limits the y-coordinate. In order to write these constraints in terms the coordinate updates, being estimated from the optimisation process, equations (4a) and (4b) are written in the following forms:

According to equations (4c) and (4d), the coordinate updates are estimated in such a way that the updated coordinates remain in the specified interval (equations 4a and 4b).

These inequality constraints for all ANs can be written in the following vector form:

where

6 Directional constraints

Directional constraints keep the updated coordinates of ANs in specific directions or azimuths. Suppose that the jth AN can only move towards a point j’th only. In this case, the following formula can be used as the equation of directional constraints

where φ j ′ j is the azimuth from j’ to j.

Writing the constraint in the form presented in equation (5a) avoids any probable singularity during the optimisation process. Since our design has four ANs, one directional constraint can be considered for each one of them. The constraints need to be written in terms of the coordinate updates; to do so these equations are linearised in the following matrix form:

where Δx is already defined in equation (3d) and

7 Optimisation model

The system of equation (4a) should be solved for the coordinate updates in a least-squares sense but subjected to the aforementioned constraints. Such an optimisation model is:

It should be stated again that considering all the constraints in equation (8) might not be possible in practice.

Today, there are different software for solving the optimisation problem (8). The theory of solving this problem is known (Bazaraa and Shetty 1976 or Grafarend and Sanso 1985). In this study, the Optimisation Toolbox of MATLAB is applied to solve equation (6).

8 Arlanda International Airport

The Arlanda international airport of Sweden is in the northern part of Stockholm. It is almost rectangular with a size of 4 km × 5 km and has three runways with an approximate size of about 500 m × 3,000 m each. Here, a local planar coordinate system is defined with an origin in the south-west of the area having the geodetic latitude φ = 59°37′10″ and the longitude λ = 17° 53′ 50″, and y-axis of the system is parallel to the western runway having an azimuth of about 10°, and the x-axis of the system is perpendicular to the y-axis and towards east. The Gaussian radius of curvature at the system origin is computed from its latitude at the surface of the WGS84 reference ellipsoid with a = 6,378,137 m, and e ² = 0.0068. A grid with a resolution of 40 m × 40 m is considered over the airport, it is rotated based on the azimuth of the western runway (extracted from Google Earth), and later their coordinates are transformed to the geodetic coordinates by:

where φ _i , λ _i are the geodetic coordinates of any point over the airport, R stands for the Gaussian radius of the curvature and the local system’s origin (cf. Jekeli 2012):

where N′ and M′ are respectively the well-known radius of the prime vertical and curvature of the local meridian at the point with the latitude φ.

Two runways of this airport are nearly elongated towards north with an azimuth of about 10°, and the other one is almost elongated from the east to the west with an azimuth of about 255°. From the satellite photo, Google Earth, we selected four ANs of M, N, O and P, in the runways of the airports, two in the most western one, and one in each one of the others. Our goal is to optimise the position of the ANs in such a way that the VC matrices of CNs, covering the area, are fitted to a criterion matrix, which is a diagonal matrix with equal diagonal elements to 1. In the first scenario, only the limiting search area constraints are applied, and in the second one, the directional constraints will be added in the optimisation process. To limit the search area around each AN, some bounds are required. We let M vary 500 m, i.e ±250 m in x-coordinates and +500 m and −2,000 m in y-coordinates. Again x-coordinate of N varies between −250 and 250 m but +1,000 and −2,000 in y; similarly, O varies between −250 and 250 m in the x-direction and +1,000 and −2,000 in y and finally the x-coordinate of the anchor P between −1,000 and 1,000 m, and its y-coordinate between −500 and 500 m. Figure 2 shows the photo of the Arlanda international airport taken from Google Earth. In the local 2D coordinate system, the initial positions of the ANs are shown with small red circles and their search areas with rectangles surrounding them. As we show in the figure, our goal is to have the ANs on the runways of the airport. In addition, as observed, the y-axis of the system is chosen parallel to the western runway for simplification. However, the choice of the coordinate system is not important as the design can be done based on the chosen system, and later, the whole network is georeferenced. Since the eastern and western runways have an azimuth of about 10°, then the y-axis of the system has the same azimuth.

Figure 2

The local 2D coordinate system, ANs and search areas of the ANs on the satellite photo, taken from Google Maps, of the Arlanda international airport in Sweden.

8.1 DOP of initial design

The initial positions of the ANs are shown in Figure 2, and as we see they are not in their optimal locations. The first step of any design of a network is to check the DOP of the network based on an initial design; Figure 3 is the map of the DOP values over the Arlanda international airport. The DOP values reach about 10 in the northern part of P, which is normal as there is no other AN to cover that part, but the other northern parts of the airport have DOP values less than 2, as those areas are covered from the east by M, N and O. Larger values than 4 are seen in the south of the airport, as all ANs are almost placed in the northern part of the airport.

Figure 3

DOP of the wireless security network based on the initial design of the ANs over the Arlanda international airport.

8.2 Optimal design without directional constraints

One possible design is to perform the optimisation process to fit the VC matrices of CNs to the criterion VC matrices, which is an identity matrix, meaning that all covariances should be fitted to zero and variances to 1. This is a rather a tough criterion and impossible for some CNs to reach such a DOP value, but it will force the optimisation process to deliver as closer as possible fit to the criterion. Therefore, it is normal to see that some CNs have higher DOP values than the criterion.

Figure 4 shows the map of the DOP of the ANs over the Arlanda International airport after optimisation considering only the search area constraints around ANs. As seen, the optimisation process pushes the ANs further away from the airport to obtain a good fit to the criterion matrix. The whole airport has a good coverage by the ANs, and large DOP values are seen more to the outside of the airport with the largest value reaching 4 in a small area in the north and south of P and O, respectively. Most of the area has lower values than 2 and the central part 1. However, the location of the ANs needs special attention from a practical point of view. They are located outside the airport after the optimisation process, e.g. P is amongst trees, and M, N, and O are around roads. Establishment of permanent ANs at M, N and O is not meaningful due to frequent uncontrollable movements of cars, increasing the risk of receiving reflected signals, multipath, and hardening the process of interference detection.

Figure 4

DOP of optimised wireless security network with search area constraints over Arlanda international airport.

8.3 Optimal design with directional constraints

To keep the ANs M and N in the western runway, it suffices to constraint them to azimuth 180° and 0°, respectively, according to the method presented before. The situation is the same for O in the eastern runway, which is almost parallel to the western runway, but slightly different for P in the runway having east–west elongation of an azimuth of about 255°. Since the y-axis of the system has azimuth of 10°, then a rotation of 245° would suffice to keep P along the northern runway. Establishment of such directional constraints, some help points are selected at a 1 m distance from each AN along the azimuth of the runways, and in fact, we anchor ANs to them so that they can vary only along the direction to these help points.

The DOP of the optimised network with the limiting search area and directional constraints are shown in Figure 5. Our study showed that optimisation solely with the directional constraint would push P and O much further away from the airport. This figure shows that these constraints could successfully keep the ANs along the runways. Nevertheless, comparing with the case where only search area around the ANs is limited, the optimal location of ANs is inside the airport, and seeing a larger DOP in the marginal areas of the airport is normal. According to Figure 5, the maximum value of DOP, about 5, is seen north of P, but the central area is well covered by DOP values less than 1.5.

Figure 5

DOP of optimised wireless security network with search area and directional constraints over Arlanda international airport.

9 Discussion

Eshagh (2022) performed optimisation of similar wireless security networks for the Landvetter international airport in Gothenburg based on TOA, AOA and TDOA. He concluded that TDOA were difficult to apply optimisation for three ANs at the airport, but in this study, we showed that the TDOA are applicable to four ANs. Also, the shape of the Arlanda airport is closer to a square than Landvetter airport with one runway. No directional constraint was considered in the optimisation of the security network of this airport, and there was not much freedom for the ANs to move to improve the DOP of the network. Applying search area constraints around the ANs is of vital importance in the optimal design process of such networks for interference localisation; otherwise, the ANs may move to places far from the airport area, but the fit to the criterion matrix will be good. By considering such constraints, the ANs can be kept inside the area but with the price of less fit to the VC matrix of criterion. Adding any additional constraints like directional constraints leads to even more misfits. However, with a proper selection of the location of the ANs and the DOP of the network, the value of DOP can reduce by about 50%, which is rather significant without adding any extra cost.

Resolution of the network was seen as significant in simulation studies done by Eshagh (2022), but in this study, no change in the coordinates of the ANs was observed for different resolutions. Here, four ANs with six observables were considered in the design, but in the synthetic test of Eshagh (2022), three ANs and three observables were applied. Finding the same design for all resolutions is a good signal because increasing the resolution costs a huge computational burden and the rate of converges of optimisation will be low. This is rather positive as by a low resolution, the same design is derived as a high-resolution network.

10 Conclusions

Designing an optimal configuration for a four-ANs GNSS interference localisation network requires the search domain constraints around each AN. Directional constraints are useful to keep the movements of the ANs in specific directions. Obviously, adding any constraint costs a lack of fit to the criterion matrix. The quadratic minimisation problem is the same as an ordinary least-squares solution with a good fit to the criterion matrix, but a lack of control over the location of the ANs. The optimal design for the Arlanda international network showed that the DOP of the network in the central part of the airport is low and in the order of the DOP defined by the criterion matrix, and the large values of DOP are mainly in the marginal areas and outside the airport. The resolution of the CNs has not any significant role in the optimal design of the configuration of the ANs, and a low-resolution grid is suitable enough for a proper optimal design. Finally, our optimal configuration design of these CNs over the Arlanda international airport showed that the maximum DOP of the network can reduce by 50%, which is significant and adds no extra cost. This means that by proper selection of the location of the ANs the chance of a successful localisation will be twice as high.