Optimization of Basepoint Configuration in Localization of Signal Interference Device

Today, almost all vehicles at the surface of the Earth are navigated by global navigation satellite systems (GNSSs). Electromagnetic signals, carrying satellite coordinates, are transmitted from satellites constantly toward all GNSS receivers mounted on vehicles. However, navigation of aircrafts needs more attention because any radio frequency interference in their navigation signal may lead to catastrophes for crew, passengers, and even people on the ground. Taking off and landing are the two most risky times of a flight and are close to an airport. Weather conditions, connection with the airport traffic control tower, and navigation signal are important factors for a successful landing or takeoff. Any intentional or unintentional interference in aircraft navigation and connection with a control tower might lead to serious problems and risk human life. Therefore, some sensors or receivers with the possibility of providing information about interference are needed. The geometric configuration of the basepoints, at which these sensors and receivers are placed, plays an important role in successful localization of an interference device. In this article, a method is developed and applied for optimization of such configuration around an area in such a way that any interference device can be localized with a higher precision.

Jamming and spoofing are two well-known types of signal interference. Jamming means transmitting a radio frequency signal into the same band as, or a band near to, the satellite navigation band of interest to prevent the signal to reach to the sensors or receivers, and spoofing is 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; see Humphreys et al. (2008) and Divis (2013). Because the power level of 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; see, for example, 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 unintentional cases such as a faulty TV amplifier, which jammed the global positioning system global positioning system (GPS) operation at a harbor in Monterey, California, for 37 days (Clynch et al. 2003). A small jammer that was used in a delivery van disrupted the ground-based augmentation system (GBAS) system aiding aircraft approaches at Newark Airport while driving on a nearby highway in 2009 (Hambling 2011; Pullen et al. 2012; Warburton and Tedeschi 2011). The Central Radio Management Office of South Korea reported several disruptions from 2010 to 2012 due to GPS jammers (Seo and Kim 2013). In Australia, Balaei et al. (2007) detected some interference and in Italy Motella et al. (2008) found some from TV signals in the GNSS band, disrupting GPS. Recognition of an interference signal among all scattered signals is a complicated process and requires skills in signal processing, which is outside the scope of this paper.

Location of an interference device can be determined from basepoints with sensors or receivers, which are able to detect and analyze an interference signal. The time of arrival (TOA), angle of arrival (AOA), and time difference of arrival (TDOA) are known measurements being used for localization. Drake and Dogancay (2004) performed the localization problem on prolate spheroidal coordinates and stated that the mathematical equations of the TDOA will be greatly simplified by them. The equations will be linear corresponding to the hyperbolic asymptotes of the TDOA. Ananthasubramanian and Madlhow (2008) investigated AOA measurements and developed a sequential algorithm; they concluded that the localization error is proportional to the AOA error variance and coverage area, and reducible by increasing the number of estimates. Thompson et al. (2009) studied the configuration of the sensors’ positions for localization of interfering devices, presented a method using differences of received signal strength measurements, and concluded that they can be alternatives to the TOA, AOA, and TDOA. Thompson (2013) investigated the interference device detection and localization by analyzing the dilution of precision (DOP), from received signal strength, and concluded that TDOA is superior to the received signal strength measurements. Numerous literature studies exist regarding localization of GNSS interference using the mentioned quantities; however, optimization of configuration of the basepoints is a novel idea, which is presented in this paper.

In geodetic network optimization, one of the purposes is to determine the optimal configuration of networks for maximizing the precision and reliability because the configuration has significant influence on the quality of the network; for more details, see Koch (1982, 1985), Xu (1989), Kuang (1996), Eshagh and Kiamehr (2007), and Eshagh and Alizadeh-Khameneh (2014). Stochastic and evolutionary methods such as genetic algorithm (GA) (Batilović et al. 2021), particle swarm optimization (PSO) (Yetkin et al. 2009; Singh et al. 2016), generalized particle swarm optimization (GPSO) (Batilović et al. 2022), and simulated annealing (SA) (Berné and Baselga 2004) are applicable for geodetic network optimization as well.

The configuration, or the geometry formed by the basepoints for localization, is also a type of geodetic network. However, in the geodetic network optimization, the configuration is optimized by varying the unknown control points in such a way that the desired precision for these points is achieved, while in optimizing a wireless network for localization of interference device the configuration of the known basepoints is optimized by varying the basepoints to reach the desired precision for the control points. So far, no optimization has been used to determine optimal configuration of wireless security networks around airports, harbors, or important infrastructures.

In this article, a new approach is proposed to determine the optimal configuration of the basepoints, at which the sensors or receivers are placed. The TOA, AOA, and TDOA of signals to these sensors and receivers are considered as observables and a criterion matrix is selected for required precision of the points over the control area. The basepoints move until the estimated variance-covariance (VC) matrices of the control points are fitted to the criterion matrix subjected to some constraints. A grid of points is designed over the control area and each one of its cells is the probable location of an interfering device. The location of basepoints is determined in such a way that the location errors of the grid points are minimized based on the type of observables and fit to the predefined criterion matrix for the grid points.

Precision of a network depends on both the quality of measurements and the network configuration. In this section, a two-dimensional control network for localization is defined and after that the observables TOA, AOA, and TDOA as well as their mathematical models are presented.

First, assume that the basepoints have known coordinates and we want to see how large the errors of the control points are over the area based on our observables. However, because the observables have nonlinear mathematical models, they need to be linearized. Generally, a linearized model is of Gauss-Markov type for a control point over the area (see Koch 2010)where $\mathbf{A}$ = coefficient matrix, whose elements are partial derivatives of observables with respect to coordinates of the control point; $\mathbf{x}$ = vector of the coordinate updates to their initial values; $\mathbf{L}$ = vector of difference between actual and computed observations from the initial coordinates; $\mathbf{\epsilon }$ = vector of random errors with $\mathrm{E}\{\mathbf{\epsilon }\}=0$, where $\mathrm{E}\{\}$ stands for the statistical expectation; ${\mathbf{C}}_{\mathbf{L}}$ = VC matrix of the observations; $\mathbf{Q}$ = cofactor matrix; and finally ${\sigma }_0^2$ = a priori variance of unit weight.

In Eq. (2a), the only parameter that carries the geometrical properties of the basepoints and the control point is the matrix $\mathbf{A}$. This matrix contains the partial derivatives of observables with respect to the coordinates of the unknown control point, which are, in fact, sensitivities of the observables with respect to the coordinates. For example, in the case of 2D localization of one interference device from $m$ basepoints, the structure of matrix $\mathbf{A}$ for TOA, AOA, or TDOA iswhere the elements of the matrix $\mathbf{A}$ are

The least-squares solution of Eq. (2a) is (e.g., Cooper 1987)and the VC matrix of the estimated coordinates

No observation is needed to compute the VC matrix of the estimated coordinates and the configuration can be optimized using Eq. (2g) and some a priori values of the measurement quality.

Eq. (2a) describes a system of equations connecting the observables to two-dimensional coordinates of an unknown control point. Obviously, at least two observables from two different basepoints are needed to have a unique solution for the system. In the case of having more basepoints, the least-squares method is applicable.

Suppose that three basepoints and three observables are used, and then the system in Eq. (2a) will have three equations and two unknowns for one control point. Our goal is to optimize the configuration of these basepoints so they cover the control area. Consequently, many such points should be considered over this area, which are probable locations of an interfering device.

As shown, the VC matrices of control points are determined from the observables and location of the basepoints. The coordinates of these control points are known as the size of the grid and its cells are designed by us. Then the basepoints are considered as unknowns but an initial configuration is needed, which can be derived from the local coordinate system defined by the grid. Therefore, initial VC matrices for all control points are computed and after that the configuration is optimized by varying the coordinates of the basepoints.

A criterion matrix is needed for the optimization process, which should be a priori known from the expected variances of the control points. In a two-dimensional network, each VC matrix has four elements, two variances as diagonal elements and two equal covariances as off-diagonal ones. Normally, a diagonal matrix with the desired variances for the coordinates of the unknown points on its diagonal elements is considered as the criterion. The initial VC matrices should be fitted to this criterion in a least-squares sense by varying the coordinates of the basepoints. In other words, by optimization the locations of the basepoints are determined in such a way that the VC matrix of the control points is fitted to the criterion.

To present this idea mathematically, the VC matrix of one point [Eq. (2g)] is expanded for the basepoints by the Taylor series (e.g., Xu 1989; Koch 1982, 1985; Kuang 1996)where $\mathbf{A}$ = initial coefficient matrix derived from the approximate positions of the basepoints; $\mathrm{\Delta }{x}_j$ and $\mathrm{\Delta }{y}_j$ = coordinate updates to the basepoints’ positions; and $m$ = number of basepoints. It will not be difficult to show that (see also Kuang 1996; Eshagh and Alizadeh-Khameneh 2014)where the structure of $(\partial \mathbf{A})/(\partial x_j)$ or $(\partial \mathbf{A})/(\partial y_j)$, which we show that by $(\partial \mathbf{A})/[\partial x_j(y_j)]$ for the TOA and AOA measurement in a 2D case of $m$ basepoints and one interference device, will be

However, when TDOA is applied this structure for the first three basepoints is

The elements of these matrices are

${\mathbf{C}}_{\widehat{\mathbf{x}}}$ on the left-hand side of Eq. (3a) is the criterion matrix, and ${\sigma }_0^2{({\mathbf{A}}^{\mathrm{T}}{\mathbf{Q}}^{-1}\mathbf{A})}^{-1}$, which is the initial VC matrix, should be fitted to ${\mathbf{C}}_{\widehat{\mathbf{x}}}$. Let us rewrite Eq. (3a) in the following form:

${\mathbf{C}}_{\widehat{\mathbf{x}}}$ is a $2\times 2$ matrix having the variances of the $x$- and $y$-coordinates on its diagonal elements and covariance between them on the off-diagonal. Eq. (3l) means that the updates $\mathrm{\Delta }{x}_j$ and $\mathrm{\Delta }{y}_j$ to the initial coordinates $x_j$ and $y_j$ are estimated in such a way that the updated ${\mathbf{C}}_{\mathbf{x}}^0$ is fitted to the desired ${\mathbf{C}}_{\widehat{\mathbf{x}}}$. Because in this study three basepoints are considered and each has two coordinates ($x$ and $y$), six unknown parameters exist in the system of equations. Our VC matrices are $2\times 2$, meaning they contain four elements and our system has four equations for each control points. However, there are many points already designed over the control area, and each one adds four equations to our systems, while the number of the unknown parameters remains constant. In this case, a huge system of equations is created and the coordinates of the basepoints are estimated in such a way that the best fit to all criterion matrices for the control points are achieved.

Let us present this system in the following form:where ${\mathbf{\epsilon }}^{\prime }$ = vector of residuals; $\mathbf{\Delta }\mathbf{L}$ = vector of differences between the elements of the criterion and initial VC matrices; $\mathbf{B}$ = coefficients matrix containing the partial derivatives of the VC matrix with respect to the basepoints’ coordinates; and $\mathbf{\Delta }\mathbf{x}$ = basepoints’ coordinate updates. The mathematical descriptions of them arewhere operator vec inserts the columns of the VC matrices below each other and converts the $2\times 2$ matrices to $4\times 1$ vectors; $({)}^{\mathrm{T}}$ stands for transposition operator of matrix algebra; and $n$ = number of control points.

At first glance, one expects to solve the system of Eq. (4a) by the least-squares method. However, unrealistic results may be obtained, for example, some basepoints may be located far from the area or they may be located close to the center and create a weak geometrical configuration. Hence, some constraints are needed for solving this system of equations, like:

The use all these constraints depends on the study area and how much freedom our basepoints must move during the optimization process. These constraints act like inner constraints (e.g., Cooper 1987; Tan 2005) that are applied in geodetic network adjustment for overcoming rank deficiency of the established system of equations due to datum defect. However, our system does not have this deficiency because the grid points already have known coordinates. These constraints are used to limit the movements of the basepoints.

To keep the basepoints outside the control area, some inequality constraints are applied for limiting the search domain of the coordinate variations; such inequality constraints are (compare Kuang 1996; Eshagh and Alizadeh-Khameneh 2014)where $w_j^{\mathrm{L}}$ and $w_j^{\mathrm{U}}$ = lower and upper bounds of the inequality that contain limiting $x$-coordinates of the $j$th basepoint, respectively; and $v_j^{\mathrm{L}}$ and $v_j^{\mathrm{U}}$ = corresponding limits for the $y$-coordinates. To write these constraints in terms of the coordinate updates being estimated from the optimization process, Eqs. (5a) and (5b) are written in the following forms (compare Eshagh and Alizadeh-Khameneh 2014):

According to Eqs. (5c) and (5d), the coordinate updates in such a way that the coordinates remain in the specified interval [Eqs. (5a) and (5b)] outside the area.

These inequality constraints can be written in the following vector form:where

To keep the average of the basepoints’ coordinates at the center of control area, the following equations can be used (compare Cooper 1987; Kuang 1996):where $x_G$ and $y_G$ = coordinates of the center of the control area, which can be determined from the control points; and $m=3$ is the number of the basepoints.

Because the nonlinear mathematical models are linearized in our optimization process, $\mathrm{\Delta }{x}_j$ and $\mathrm{\Delta }{y}_j$ are estimated and not the coordinates directly. Therefore, Eqs. (6a) and (6b) need to be rewritten in terms of these coordinate updatesorwhere

Bearing angles from the center of the area to the basepoints can be used for making a constraint for distributing the basepoints around the area. By assuming that the summation of all these bearing angles is zero, we can write (compare Kuang 1996)

Again, Eq. (7a) should be written in terms of the coordinate updates; therefore, its linearized form will be appliedwhere $l_{jG}$ = distance from the center of control area and the basepoint $j$. Eq. (7b) can be written in the following matrix form (compare Kuang 1996):where

The system of equations in Eq. (4a) should be solved for the coordinate updates in a least-squares sense but subjected to the mentioned constraints. Such an optimization model is (compare Koch 1982, 1985; Kuang 1996)

Considering all these constraints, Eq. (8) might not be possible in practice. As shown in the synthetic test, all basepoints have the freedom to move, but their movements can be limited by the three constraints so they can move around the area using the search area constraint and around the center of the area with a good distribution. However, in the real study for the international Landvetter Airport, based on the presence of forests around the area, keeping the basepoints around the area is not practically possible and only the limiting search area around each basepoint is suitable to keep them in area.

Today, there are different software programs for solving the optimization problem in Eq. (8). The theory of solving this problem is known (e.g., Bazaraa and Shetty 1976; Grafarend and Sanso 1985). However, the optimization problem in Eq. (8) is not global because of the involvement of different constraints. The unconstrained optimization may diverge, or the basepoints may move outside the area or toward themselves or become colinear. Applying at least the search area constraint is a necessity for obtaining a convergent solution. The resolution of the grid of control points plays a role in the rate of convergence and computational time. In addition, a threshold should be a prior defined to stop the iterative optimization; in this study, the norm of the coordinates updates is computed and when it becomes smaller than 10 cm during the iterations, the optimization is stopped.

To test the presented theory, a rectangular area of $4\times 5\mathrm{km}$ is considered with resolutions of 10, 20, 40, and 80 m. Three basepoints, M, N, and O, are chosen manually around the area; see their initial coordinates in Table 1 and their locations in Fig. 2(a) by small circles. Fig. 2(a) shows the pattern of the square root of the trace of the VC matrix of the control points, i.e., horizontal dilution of precision (HDOP), based on the TOA observables and a resolution of $40\times 40\mathrm{m}$. Because the HDOP patterns for different resolutions were similar, only one is presented. Generally, the HDOP based of the TOA has no units, but we considered a priori precision of 1 m for the TOA distances to give a unit to the computed HDOP. Such a precision seems to be large for geodetic networks, but it is acceptable in localization of the interference device by wireless networks (e.g., Trinkle et al. 2012). As observed, the largest error reaches to about 2.2 m, close the farthest point O at the upper-right corner of the area. A diagonal criterion matrix with equal diagonal elements to 2 m was considered for fitting the VC matrices of the control points to it. The basepoints were kept 400 m, as an example, outside the area from the margins by using the inequality constraint in Eq. (5e). Fig. 2(b) shows the optimized location of the basepoints on the map of the HDOP. A significant reduction of HDOP is seen with the maximum reducing to 1.6 m. In addition, the configuration forms almost an equilateral triangle, which is more stable than the initial design forming an isosceles triangle. However, the weak points having larger HDOPs are seen around the basepoints and elongated toward the center of the area. This could be expected because two of the three distances from each control point and to the three basepoints are considerably longer than the other one, forming a weaker geometry than the points located at the center of the area.

A similar HDOP pattern is seen when the AOA are considered. The same criterion matrix was considered for the control points and a priori variance factor of 60 arc-seconds as the accuracy the measured AOAs. Fig. 2(c) shows that the HDOP before optimization reaches to about 3.5 m, and after optimization, Fig. 2(d), it significantly reduced to about 1.8 m with a symmetric configuration.

Fig. 2(e) is a similar HDOP for the TDOA, with a priori variance factor of 1 m, based on the same resolution before optimization. It shows that the HDOP reaches to about 3.5 m and the weak points are seen at the corners. Generally, a better coverage over the area can be seen using the TDOA. The HDOP after optimization is shown in Fig. 2(f). As observed, a symmetric geometry is created by the optimization process and the weak points are around the corners of the area, but with a proper coverage over the inner area meanwhile reducing the HDOP to 2.5 m.

The HDOP derived from the TDOA gives a good coverage over the area, and the weak points have large HDOPs and are not close to the basepoints. However, the mathematical models are more complicated because of their differential nature.

To test sensitivity of our optimization method based on the observables, the process was repeated with different resolutions. Table 1 gives the optimal horizontal positions of the three basepoints after optimization, while the initial positions of these points were M (0.5, 2,500.5), N (2,000.5, 0.5), and O (4,000.5, 5,000.5) prior to optimization.

According to Table 1, when the TOA is used, the $y$-coordinate of M remains the same for resolutions $10\times 10\mathrm{m}$ and $20\times 20\mathrm{m}$, but the $x$-coordinate remains constant for all. There are variations in the $x$-coordinate of N but it is almost constant for these resolutions as is the $x$-coordinate of O. Generally, with resolutions less than $20\times 20\mathrm{m}$, the optimized coordinates of the points are almost the same as the observable.

Similar results are obtained when the AOA is used. A resolution of $30\times 30\mathrm{m}$ is good enough for optimization, but for TDOA no convergence is seen by increasing the resolution. Nevertheless, the differences between the coordinates obtained from optimization using the two successive resolutions are in order of the resolutions. Therefore, the TDOA is more sensitive to the selected resolution and no specific resolution can be determined for optimization based on it.

Generally, applying a dense grid for the control points is suggested, but it takes a huge computational burden and long time to perform the optimization process. One important issue that needs to be highlighted is that our study area is rectangular $4\times 5\mathrm{km}$, close to a square.

So far, our simulation study showed that the presented theory works with some constraints. Designing a security network is not as simple as what was done in the simulation study unless the area is flat with no signal hindrance. Here, an international airport is considered and positions of three basepoints are optimized in the designed control network with the assumption of the presence of the interference device inside the airport. Not all our proposed constraints can be applied for this study area. First, there are forests around, which are not suitable places for our sensors. Therefore, spreading the basepoint outside and around the airport is not realistic. Second, search areas around the basepoints are limited because there is not much freedom for them to move.

The international Landvetter Airport of Sweden in Gothenburg was selected, and a local planar coordinate system was defined for the airport. The Gaussian radius of curvature at the system origin was computed from its latitude at the surface of the WGS84 reference ellipsoid with $a=\mathrm{6,378,137}\mathrm{m}$, and $e^2=0.0068$. A grid with a resolution of $40\times 40\mathrm{m}$ was considered over the airport; was rotated based on azimuth of the most left takeoff way (extracted from Google Maps), and later their coordinates were transformed to the geodetic coordinateswhere ${\phi }_i$ and ${\lambda }_i$ = geodetic coordinates of the $i$th point of the grid; ${\phi }_O$ and ${\lambda }_O$ = geodetic coordinates of the local system origin; $x_i$ and $y_i$ = local Cartesian coordinates of the point; and $R$ = Gaussian radius of the curvature and the local system’s origin (compare Vanicek and Krakiwsky 1986; Jekeli 2012)where $N$ and $M$ = well-known radius of the prime vertical and curvature of the local meridian at the point with the latitude ${\phi }_O$, respectively.

Landvetter Airport’s runway is elongated in southwest–northeast direction with an approximate bearing angle of 26°. The area is almost rectangular, with length of 5 km and width of 2 km. Based on the geometry of the airport, no symmetric configuration can be considered for the basepoints. Therefore, from the satellite photo and Google Earth, we selected three suitable places; Fig. 3 shows thee selected areas by the squares for optimal locations of the basepoint. In fact, an initial position was considered for each basepoint, which varied in a square with the size of 400 m.

Fig. 4(a) shows the maps of HDOP for the TOA over Landvetter Airport and the small circles represent the basepoints before optimization. The points having larger HDOPs were extended along the runway of the airport. The HDOP reaches to about 6 m over the area. Fig. 4(b) illustrates the optimized positions for the basepoints, and it shows that the optimization has reduced the HDOP to 4 m with slight change in the initial positions; see Table 2. Fig. 4(c) is a similar map to Fig. 4(a) but for the AOA with the same initial positions for the basepoints. After optimization, the map in Fig. 4(d) is obtained, and small changes in the positions of the basepoints lead to reduction of the HDOP to less than 4 m. Because of the special narrow rectangular form of the airport, the initial design based on the TDOA was not successful; the reason could be that the hyperbolic curves, which are created from each two points, do not make a good geometry. The HDOP at some points became extremely large and unrealistic, and even if the optimization was performed the constraints were violated and the basepoints moved outside the area. For this reason, no result is presented for design using the TDOA.

Table 2 lists the initial and optimized positions of the basepoints M, N, and O. Optimization based on the TOA and AOA leads to the same positions for the basepoints M and O, but N is more influenced by the type of observables. Generally, M and O are closer to each other, but N is far from both. These points form a triangle having a small angle at N. It is known that an equilateral or equiangular triangle is a strong geometric form. When a triangle has a different form, it has weaker geometry, such as when M, N, and O form a triangle having such a small angle at N. This is the reason N moves during the optimization process, because it has direct influence on the strength of the geometry of the basepoints.

TOA, AOA, or TDOA of signals into sensors or receivers can be used for localization of a signal interference device. This article showed that the optimal geometrical configuration formed by basepoints, at which these signal arrivals are measured, strongly depends on the type of arrival and the ability of sensors or receivers to measure these arrivals. A design using the AOA showed a better stability and robustness with respect to the resolution of the designed grid compared to the rest of its rivals, and optimization based on the TDOA improves coverage over the control area. Generally, our optimization method needs some constraints for stabilizing the solution; otherwise, the basepoints may form improper configuration during the iterative optimization process, e.g., they may move toward each other or become collinear. However, in our real case study, the search domains for the basepoints were limited and subjecting the optimization process to extra constraints, as done in the simulation test, led to infeasible solutions. In other words, the choice of the constraints is highly dependent on the shape of the area.

The data that support the findings of this study are available from the author upon reasonable request.

The author is thankful to both reviewers and their fruitful comments enhancing the quality of the paper.