Passive Radar Transmitter Localization Using a Planar Approximation

In this article, we review and compare different methods for the localization of an illuminator of opportunity for passive radar systems employing cooperative targets. For assessment and comparison, we use numerical simulations and real-world radar measurements. To achieve reliable operation under real-world conditions, a correct and robust association of the bistatic radar target observations with ground truth information is necessary. We introduce a novel solution to this issue with an efficient preprocessing algorithm based on a planar approximation.


I. INTRODUCTION
T HIS work addresses the problem of finding illuminators of opportunity for a versatile, mobile passive radar system. In such a bistatic system, information on the transmitter positions is essential for the interpretation of the measured range-Doppler data. In most passive radar systems described in the literature, the positions of these noncooperative 1 illuminations are assumed to be known [1]- [3]. Various databases exist, from which such information can be retrieved, but required data such as location, height, transmitter frequencies, and power may have to be compiled from different sources. Such public databases are often organized as private projects and, therefore, do not guarantee the actuality, accuracy, or completeness of the data they contain. Furthermore, commercial use is often not permitted. It would, therefore, be desirable to use a truly independent system that does not rely on such external information. We present a method for the localization of such suitable (see [4]) illuminators, utilizing cooperative targets that transmit their ground truth information. Such a procedure was first described by Yi et al. [5].
Malanowski et al. [6], show one possible implementation of such a technique. In this article, we present alternative approaches to do so and examine their performance with numerical simulations, especially with respect to the standard deviation of the measurement error, which Malanowski identified as the main issue in the loss of performance when comparing real-world results to the simulation results of Yi et al. [5].
In real-world scenarios, we usually have to deal with multitarget situations. Therefore, a correct solution (or minimization) of the equations for the illuminator position requires that the observed range and Doppler information (from the radar measurements) are correctly assigned to the respective cooperative targets (with their known positions and velocities). Observations may even be discarded if the correct target cannot be assigned or the target ground truth is not known. Due to the complicated dependencies of bistatic range and Doppler from the 3-D positions of the receiver, target, and (unknown) illuminator, and from the target velocity vector, this assignment has no obvious or direct solution. With the planar approximation (PA) method, we show a novel simplified approach for the association of the bistatic measurement results to the ground truth information, which is computationally efficient and which would be suitable for real-time implementation.

II. SYSTEM DESCRIPTION
The measurements described in this article were conducted with the ALFA passive detection system, which is a software-defined radio (SDR) based multiband receiver system [7]. The system is equipped with a circular eight-element VHF array, 2 among others, which is used for Frequency modulation (FM)-radio-based passive radar. The basic principle of such a system is depicted in Fig. 1; the setup consists of a receiver B and an illuminator I. Similar to a conventional bistatic radar, it produces range-Doppler measurements for the target T . The bistatic range R is calculated from the time difference of arrival between the reflected signal and the line-of-sight signal tap. It is defined as The Doppler frequency is calculated from the frequency shift resulting from the target's velocity v and the center frequency of the illuminating signal f I . The Doppler frequency D for this scenario consists of the dot product of v with the difference of two radial components, one from the Illuminator I to the target T and one from the target to the observer B. It is defined by Our implemented system does not employ a dedicated reference antenna for the direct signal tap; instead, a delayand-sum beamformer [9, pp. 23-32] is used to form multiple virtual directional antennas. One virtual antenna points toward the transmitter of opportunity, which is used to generate the reference signal. In total, 16 additional virtual antennas point in 16 equidistant directions, covering the full 360 • in azimuth. Fig. 2 depicts the signal processing for the passive radar application. First multipath effects (i.e., clutter) are removed from the reference signal, using the constant modulus algorithm [10]- [12]. In order to reduce the direct signal component in the 16 surveillance channels, the normalized least mean squares algorithm is used to subtract the cleaned reference signal from them [13], [14]. The surveillance channels will then primarily contain reflections. Each surveillance channel signal is split into smaller pieces, which are individually correlated with the reference channel and filled into a matrix. An fast Fourier transform (FFT) over the slow-time axis of the matrix yields the range-Doppler matrix. A constant false alarm rate (CFAR) algorithm is used to detect the targets in the range-Doppler matrix [15]. The resolution of the detected targets is enhanced by reconstruction of the signal in the range and Doppler domain using the sinc-interpolation method [16], which is applied to all 16 virtual antennas. As a result, each target is represented by 16 amplitude values versus observation angle of the respective virtual antenna. To enhance the angular accuracy, this signal is in turn reconstructed with sinc-interpolation and its maximum is evaluated as the angle toward the target. 3 With this process, a list of targets are created for each illuminator.

III. GROUND TRUTH ACQUISITION
For our work, we use automatic dependent surveillancebroadcast (ADS-B) to acquire the ground truth information of the cooperative targets [19]. Since ADS-B transmitters are installed on most commercial aircraft, we found a multitude of cooperative targets at all our test sites so far. The ADS-B signals are received with an additional RTL2832U based SDR independently from the passive radar system, and all incoming messages are logged while passive radar measurements are taken. It has to be noted that the ADS-B messages are not synchronized to the acquired passive radar measurements, therefore, the ground truth information cannot be mapped to the range-Doppler targets directly. For the mapping, the position and speed vectors of the targets are interpolated with the two nearest available ADS-B positions. We use a simple linear interpolation, which assumes that the aircraft is not maneuvering between two ADS-B updates. Of course, the unknown error resulting from this assumption adds further to the range-Doppler measurement error.

IV. IMPLEMENTATION OF THE ILLUMINATOR FINDER
We assume a dataset for each illuminator containing range-Doppler measurements with the corresponding ground truth from multiple cooperative targets, with several data points each. A procedure for assigning the range-Doppler measurements to the ground truth is described later. By applying a dataset, with a total of N measurements, to (1) or (2), we can describe the cost functions for range and Doppler Figs. 3 and 4 show a plot of these functions for an area of 200 km by 200 km around the receiver. For both plots, the same set of randomly generated flights (black dots) are used. It is obvious that both the Doppler data and the range data 4 have a clear minimum at the position of the illuminator. Solving these equations for the illuminator position leads to an overdetermined nonlinear system of equations, for which we will present three different solution approaches in further detail.

A. Gauss-Newton Approximation
One possibility for solving the equations for the illuminator position is the Gauss-Newton iteration (GNI) algorithm [20, p. 254 ff.], which is an iterative procedure for solving nonlinear least-squares (LSs) problems [20, p. 245 ff.]. The general iteration of the algorithm is described by where r describes the residuum between the measured value of either range or Doppler and their calculated values. J is the corresponding Jacobian matrix, 5 and γ is a step width (0 < γ < 1) to ensure the convergence of the algorithm. This type of iterative solution furthermore needs an initial guess of the solution as a starting point. We show two implementations of this method, one for range values and one for Doppler. The first step in both is to move the receiver coordinates to the origin of the coordinate system by subtracting B from T . We get with I = (I x , I y , I z ) T and O = T − B, to describe the range of a target. The iteration is described by The residuum is the difference between the calculated range value and the measured range value, which leads to the residuum vector J is given by The equations for the Doppler method can be developed in a similar way, starting from To evaluate the reliability of this approach and to find reasonable values for γ , this approach is tested with simulated data. 6 Furthermore, the influence of the number of targets on the reliability is examined. For each combination of γ and measurement count, 10 000 simulations are evaluated. A simulation run is counted as successful if the GNI reaches the illuminator position with a distance smaller than 1 m. Otherwise, the simulation is counted as failed. Fig. 5 shows the percentage of failed simulations. It can be seen that the error quota does not depend on the number of targets used in the simulation and that a γ value below 0.5 yields the best results. The data furthermore show that even with ideal, simulated data without added error, the algorithm fails in a small percentage of cases to find the correct illuminator position. This behavior can be described by local minima in nonlinear optimization.

B. Spherical-Interpolation (SI) Method
Another method to determine the illuminator position by using (6) is shown in [21] and is called the SI method. The difference between Malanowski's algorithm in [21] and 6 For the simulation, the receiver and the illuminator are defined at fixed positions, 10 km apart. The position of the targets and the initial guess for the illuminator position are chosen at random for each simulation. The range and Doppler values are calculated with (1) and (2). No additional measurement error is added to the calculated values. The target position T = (T x , T y , T z ) T is chosen with random cartesian coordinates around the receiver position B = (0, 0, 0) T , with T x , T y and T z each between −10 km and 10 km. our approach is that the target and the illuminator have interchanged their roles. For a single measurement n, the target position is described by the vector The corresponding range is denoted as R n . Before the SI method can be applied, (6) has to be rearranged to with w n = R n − O 2 x,n + O 2 y,n + O 2 z,n , the unknown position of the illuminator I = (I x , I y , I z ) T and For a dataset with a total number of N measurements, (9) can be written in matrix vector notation as For N > 3, (11) describes an overdetermined system of equations. We now assume that the value of P is known and define an error vector which has to be minimized with respect to I in terms of LS. The solution to this minimization problem is not exact but unambiguous and yields [21] with the pseudoinverse matrix S + . To get a solution for P, we substitute (13) into (12), which results in a new error vector [21] P = (I − SS + )( V + 2P W ) which now has to be minimized with respect to P. The solution forP is given by [21] with T = I − SS + . Now, (6) can be solved by calculating P with (14) and using this result in (13) to get the illuminator position˜ I. It is clear that, compared to the GNI approximation (see Section IV-A), this represents a closed solution to the problem, so no initial guess is necessary for this procedure.

C. Spherical-Intersection (SX) Method
In [21], Malanowski describes the SX method, another variant of a closed-form solution. Starting from (13), we substitute If we now substitute I T I = || I|| 2 and insert this into (10), we get the quadratic equation The two solutions of (16) are given by With (17) and (13), we can calculate two possible locations for the illuminator. The correct solution can be chosen by comparison of the value of the cost function (3); the position with the lowest cost is used [21]. If due to larger measurement errors, the discriminant of (16) is less than zero, the solution of (17) has an imaginary component. In this case, only the real component of the solution is used [21].

D. Initial Dataset Matching Using PA
The approaches described earlier require a dataset with range-Doppler measurements associated with a ground truth coordinate. A possible solution to solve this target association ambiguity is proposed by Malanowski in [6]; the proposed method checks for each combination of bistatic measurements with available ground truth information the Fig. 6. Overview of the used coordinate systems. The ECEF coordinate system is marked by X , Y , Z. The local system, tangential to the surface of the earth is marked as X ', Y ', Z'.
minimum value of the cost function. However, for larger datasets, this procedure becomes extremely computationally complex and is, therefore, no longer applicable. To solve this problem, we are suggesting a new approach.
In our system, the matching of a range-Doppler target to the ground truth coordinates is done by comparison of the true bearing, calculated from the target coordinate, and the bearing extracted from the beamforming process. To match a combination, we allow for a deviation of ± 2.5 • . This simple approach manages to filter out a lot of false combinations, but some unwanted, incorrect assignments are still left. These false associations lead to problems with the calculation of the illuminator position, which are described in Section VI.
What we need at this point is a good approximation of the illuminator position to filter out the incorrectly assigned measurement results from our dataset. To this end, we conceived the PA method as a robust procedure simplifying the problem to a single-variable problem by further utilizing the measured bearings toward the illuminators. 7 For that purpose, we assume that the illuminator is on a tangent plane to the earth originating at the receiver location, thereby ignoring its actual height above ground and the curvature of the earth. We can describe this plane with a coordinate system, whose z-axis is the extension of the vector B, as can be seen in Fig. 6. The y-axis of this system points toward the North Pole. In this coordinate system, the position of the illuminator I can be described with the measured bearing α and its distance R With this simple approach, we can calculate an illuminator position for every single measurement we have, without the necessity for combining multiple measurements, which may bear the risk of adding false assignments to our dataset. Incorrectly assigned range-Doppler measurements would lead to outliers of the calculated position on the straight line given by the angle toward the illuminator. The correctly assigned measurement results appear distributed around the actual location of the illuminator. To minimize the outliers' influence, the transmitter position can now be formed from this distribution with the median instead of simply calculating the mean value.
Since the position of the targets and the receiver location is described in a global earth-centered earth-fixed (ECEF) coordinate system, we decided to transform this local coordinate system to the global one. The transformation is described by with the unit vectors The angles θ and ϕ are described by with arctan2 being defined in the usual way. 8 Applying (20) to (19) results in By applying (22) to (6), we get an equation for the range with only one unknown variable, R. This equation R(R) can be solved using (5) to calculate a value R for each measurement, resulting in an illuminator position using (22). By taking the median of all these calculated positions, quite a good guess of the illuminator's position can be made. Fig. 7 shows an example of this process, using a normally distributed range error with a standard deviation of σ = 300 m. First, this guess of the illuminator's position is used to filter the dataset further and remove false combinations of range-Doppler and ground truth measurements. This is done using (6) and (8) to calculate the range and Doppler values, using the estimated illuminator position and the ground truth data. The results are compared to the measured range and Doppler values from the radar. If they differ too much, the combination is removed. Second, the illuminator position is used as the initial guess for the GNI approximation, which now uses the double-filtered dataset. As a result of this more precisely known starting position, the number of converging solutions increases.

V. EVALUATION OF THE METHODS WITH SIMU-LATED DATA
In this section, we will evaluate the different presented methods based on randomly simulated scenarios. 9 First, we examine the effect of using the position calculated with the PA as an initial guess for the GNI-based methods. For that purpose, 10 000 simulations for each of the four combinations, GNI either based on range or Doppler and using the initial guess from the PA or a random position, are conducted. To take measurement error in the angular estimation into account, a uniformly distributed random error of ±2 • is added to the bearing α. No additional error is added to the range-Doppler data for this examination. Table I shows the results. 10 The results shown make it clear that, although the error probability is low, the standalone GNI approaches are unreliable as the errors can exceed 9 For all simulations, the illuminator position is assumed at a fixed position, 35 km away from the receiver with the antenna position 100 m above the ground. The angle α describes the bearing from the receiver toward the illuminator. The targets are simulated as commercial airliners, with linear flight path intersecting with the measurement area. To simulate these flights, the start and end positions are chosen at random, both in a range between 50 km and 80 km around the receiver, with an altitude between 2 km and 12 km. The flights are modeled as a linear movement with constant velocity from the start point to the endpoint. The velocity is chosen randomly between 214 m/s and 286 m/s, resembling typical airliner cruising speeds. To get viable distances between the measurements, a time step of 30 s is used. For each simulated dataset, flights are generated until a set of 100 measurements are complete. A single measurement consists of the target position, its speed, and the range-Doppler value. 10 The error already described in section IV-A appears in the error quota, which is again defined as the percentage of calculations with an error bigger than 1 m.  10 km. These errors are mitigated completely by combining the GNI method with PA, thus severely improving their reliability.
In the next step, we examine the prior neglected influence of the range/Doppler measurement error. Therefore, a normal-distributed error is added to the calculated range and Doppler values from the simulated radar. The standard deviations σ R and σ D of this error are varied in this investigation, and for each error distribution, 1000 simulations are conducted. The added error of the angular measurement is still ±2 • . The distance from the calculated position to the actual illuminator position is evaluated as a measure of the method's performance. For better illustration, the results are shown both in histograms with two exemplary results in Figs. 8-12 and as a graph with the mean error and the standard deviation of the results in Fig. 13.
To make the results comparable, each set of simulated targets is tested with all described methods. The first method applied to the dataset is the PA. Figs. 8 and 13 show that the PA method is fairly resistant against added range error. For range errors with an error distribution with σ R > 500 m, the mean error is approximately constant. For smaller values, it is even slightly larger, which could be a result of the nonlinear median filtering. The PA results are then used as an initial guess for the GNI methods; the outcome of these investigations can be found in Figs. 9 and 10. In both cases, the mean error for the calculated illuminator position grows linearly with the added measurement error (see Fig. 13). A remarkable difference between the two Fig. 9. Error distribution for the range-based GNI with different error distributions for the measured range, with the initial guess calculated by the PA method. methods can be found in the progression of the standard deviation compared to the mean value. For the range-based methods, the standard deviation is always lower than the mean error, but for the Doppler-based methods, it is partly larger and overall closer to the value of the mean error. This means that for comparable mean errors, the results of the Doppler-based methods have a greater spread, making the range-based methods more reliable. The results for the SI and SX methods are shown in Figs. 11 and 12. The mean error increases linearly with the error on the input data. However, in both cases, the slope is considerably larger than with the GNI method. For the SI method, it is nearly doubled, for the SX method nearly tripled. The behavior of the deviation is comparable to the GNI method.  The advantages of the GNI-based methods can be explained by approximations in the closed solutions of the SI and SX methods [22], [23]. Finally, the influence of the angular error, previously modeled at ±2 • , on the PA and the dependent GNI procedures was investigated. For this purpose, simulations were performed with a range error of σ R = 480 m and a Doppler error of σ D = 2.8 Hz, respectively, and a fixed angular error was added to the angle toward the illuminator. For each variation of the angular error, 1000 simulations were conducted. In Fig. 14, it can be seen that with increasing angular error the error of the PA method also increases linearly. However, the performance of the GNI methods is not affected since the result of the PA is merely used as an initial value of those methods.
Summing up these results, the GNI, in combination with the PA method, yields the best results for all considered measurement error distributions.

VI. REAL-WORLD EXAMPLE
The methods described earlier are now tested with real-world measurements, which were carried out in The Netherlands near The Hague. An illuminator 11 transmitting FM Audio broadcast at 95.2 MHz is chosen. 12 As described 11 The respective transmitter is located at 52 • 08 13.6 N 4 • 38 47.3 E, approximately 90 m above ground. It is part of the "Zendmast Bezuidenhout" broadcast tower [24]. 12    The dataset is generated from the range-Doppler data by applying a CFAR with a 15 dB SNR threshold. The targets are then matched to interpolated ADS-B targets by the measured angle, allowing for a maximum deviation of ± 2.5 • . 13 In this case, an illuminator whose position was known was chosen to validate the results: for the sake of the experiment, it was assumed to be unknown. First, we apply the abovementioned methods directly to this dataset, which yields the results described in Table II. It is apparent that all methods, except the PA, fail to get near the illuminator position due to the abovementioned incorrectly assigned targets. Looking at the presentation of the cost functions (3) and (4) shown in Figs. 15 and 16, it becomes clear why these methods fail. For neither of the two datasets, the cost function has a minimum at the actual position of the illuminator; on the contrary, the Doppler function seems to have its maximum close to the position of the illuminator.
To filter these false associations from the dataset, we use the result of the PA method. Based on these results, we 13 The threshold value for the angle was determined empirically based on the available measurement data and probably overestimates the system's actual measurement accuracy. It is deliberately chosen small in order to preselect the dataset for the following processing strongly.   14 The results achieved with these filtered datasets are shown in Table III. The position found with the PA is again used as an initial guess for both GNI methods.
The range-data-based methods improve drastically when based on the filtered dataset, while the improvement of the Doppler-data-based method is not as impressive, as the error is comparatively large. All range-data-based methods produce an error smaller than 1 km. The best result for the real-world data is achieved by the Malanowski's SI method [21], reaching the illuminator position with an error of less than 500 m. To illustrate the improvement of the dataset, we look again at the representation of the cost function in Figs. 17 and 18. For the range data in Fig. 17, it is now clear that the minimum is at the illuminator position. For the Doppler data in Fig. 18, the position of the minimum is not so clear, but at least the illuminator position is now in a range of minimum values, even if these extend over a large range. We assume that the measurement error in the Doppler data is much greater compared to the range data. Since the error in the data is a combined error between the actual measurement error and the error within the interpolated  ADS-B data, it is not possible to analyze the origin of this error in further detail.

VII. CONCLUSION
We have presented different methods to find the position of an unknown illuminator of opportunity in a passive radar scenario using the ADS-B information of cooperative targets and evaluated their performance with simulated datasets. All of these approaches fail in a real-world scenario due to false-positive combinations of radar observations with ground truth data. To solve this problem, we devised the PA method, which yields a robust approximation of the illuminator position. We have demonstrated that this approximation can be used to prefilter a dataset containing falsely matched targets and, furthermore, can be used as an initial guess for iteration-based solving approaches, improving their reliability. We have shown that our technique improves the performance of iteration-based solving methods with simulated datasets and that it is a robust approach for solving the target association ambiguity problem in a real-world dataset.