Three-Dimensional Localization of Active Aerial Targets Using a Single Terrestrial Receiver Site

This paper proposes a method for the three-dimensional localization of an active aerial target by a single ground based sensor. The proposed method employs the time and frequency differences of arrival of the signal received directly from the aerial target and the signals received after being reflected from some large auxiliary terrestrial targets (pseudo-sensors) with known positions on the ground. Due to the terrestrial nature of the main and the pseudo sensors, it is impossible to solve for the target's altitude using traditional methods. The proposed method employs target motion analysis to obtain target position including its altitude with acceptable accuracy and low computational complexity. Presented simulations confirm acceptable accuracy of the proposed method in determining three dimensional position of the target despite limited number of the pseudo sensors and its low computational complexity.


I. INTRODUCTION
HE problem of radiation source localization has drawn considerable attention in different applications such as radar [1], sonar [2], wireless sensor networks [3], and navigation systems [4]. Time difference of arrival (TDOA) and frequency difference of arrival (FDOA) are amongst the most frequently used measurements to solve for signal source location [5]. Conventional single site passive radars use the target's radiation to determine the signal direction and its parameters. In single-site localization, in addition to the angle of arrival of the signal, the target's distance is also required. The major problem of the conventional single-site passive radars is that they cannot obtain the target's distance. To solve this problem, reflection of target signal, from one or more large targets (pseudo-sensors or auxiliary targets) with known position(s), has been proposed for source localization [6]. In passive multisensor localization scenarios, telecommunication links are used for communication between and synchronization of central and remote stations. Using telecommunication links, however, can violate covert operation of these systems. To the contrary, in single-site localization systems, no telecommunication links are used, and signal reflection from large objects is used to simulate the role of extra sensors for target localization [7]. In these methods, the cross-correlation between the reference and the reflected signals is used to extract the needed time and frequency differences for target localization [8].
In ground based single-site localization scenarios, due to the terrestrial nature of the stations, obtaining a three-dimensional position for the target is a complex and challenging task. To address this issue, the three-dimensional localization problem is transformed to a two-dimensional problem with specified hypothetical altitude, simplifying the problem to large extent, and a closed-form expression for the target's coordinates and its two-dimensional velocity are also independently obtained. The two-dimensional coordinates of the target are derived in an explicit form using a geometric method. Two of the most frequently used geometric methods for signal source localization is the lines of position (LOP) [9,10], and the location on the conic axis (LOCA) methods [11]. The use of the LOP and LOCA methods has many disadvantages including computational complexity, difficulty of error performance analysis, and inability to solve for parameters such as target speed. The geometric method utilized in this paper is a kind of analytic solution introduced by Fang [12]. Fang's method provides an exact solution, however, it does not make use of redundant measurements available by additional sensors to improve localization accuracy. To solve for the target position in an n-dimensional space, we need to have at least + 1 sensors, but these sensors should not be positioned in a sub-space with fewer than dimensions [13]. Here, however, all of the sensors used for localization are terrestrial, making this method ambiguous in determining the target altitude in three-dimensional space. In this paper, we use target motion analysis [14] to resolve this ambiguity, an idea inspired by the interacting multiple model (IMM) estimation method [15,16]. In the proposed method, probabilities are assigned to each of a number of hypothetical altitudes based on a specific criterion and then the previous altitudes are updated based on the allocated probabilities. The target's two-Three-Dimensional Localization of Active Aerial Targets Using a Single Terrestrial Receiver Site   Saber Kaviani and Fereidoon Behnia   T   > REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <   2 dimensional coordinates and velocity are obtained independently in the next iteration based on the updated altitudes. This iteration continues until an acceptable level of accuracy is achieved. By extending the knowledge of existing works, the present study contributes to resolving the ambiguity of altitude in three-dimensional localization. In addition,  We achieve a closed and explicit form for magnitude and direction of the velocity vector by assuming that the altitude of the target is known.  We present an algorithm for obtaining the correct altitude of the target. This algorithm works as follows: first, a probability is assigned to each of a number of hypothetical altitudes, and then the following locations of the target are estimated using the Kalman and an unbiased finite impulse response (UFIR) filter. This structure operates in parallel, and the weighted estimation errors, which are calculated from the difference between the predicted position in the previous time step and the target's current position, allow optimal output selection between models attributed to the target.  Finally, by presenting a specific criterion based on the target's current and predicted positions, we reallocate each altitude's probability and then determine the new altitude for the next iteration accordingly. The iterations continue until a desired level of accuracy is achieved. The organization for the rest of this paper is as follows. In Section II, we find the target's position and velocity in two dimensions by assuming that the target's altitude is known. Section III presents the proposed algorithm for determining the unambiguous altitude in three-dimensions. Simulation results and subsequent discussions are given in Section IV, and finally, conclusions and a brief summary of the advantages and disadvantages of the proposed method are presented in Section V.

II. CLOSED-FORM SOLUTION FOR THE THREE-DIMENSIONAL LOCALIZATION AT A HYPOTHETICAL ALTITUDE
The localization method presented in this study is based on measurement of time and frequency difference of arrival of the signals collected from the central terrestrial station's two receiving channels. As shown in Fig. 1, we use the echo of the source signal reflected from the known targets (pseudosensors) and the source direct path signal to calculate TDOA and FDOA. In this step, the altitude of the target is assumed to be known.
In this study, only two dominant auxiliary targets with known positions are used to solve for the target position. This is because, the method used in this paper does not make use of redundant measurements provided by additional pseudo sensors and so there is no benefit in resorting to more auxiliary targets.  Fig. 1 provides an overview of the receiver structure. The signals reflected from the auxiliary targets are cross correlated with the direct path signal, and the correct Doppler and delay values are extracted using the constant false alarm rate (CFAR) block. The reason for using the CFAR block is to provide some degree of immunity to unwanted/uncorrelated high power signals. In addition, by choosing the suitable type and setting proper coefficients for the CFAR, the effects of unequal gain of the receiver in the frequency domain can be partially compensated. Fig.2 illustrates how to use targets with known positions in determining the target's position.

Known Position
Known Position In this method, three sensors and two independent measurements are sufficient for two-dimensional positioning, and the altitude of the target is assumed to be known. After obtaining the two-dimensional coordinates of the target, we convert these coordinates into the spherical distance, the azimuth angle, and the elevation angle as below 2 2 2 In the next step, the target speed is obtained explicitly in two dimensions in a closed-form. To do end, first of all, we examine the vector form of FDOA equation and then derive its scalar and nonlinear equivalent, and then, we obtain the target position in three dimensions for a predetermined altitude. The linear and vector-form of the equation is 0 0 where is the position of the i-th sensor, indicates the position of the target, is the value of the time difference of arrival, ̇ is the value of the frequency difference of arrival and 0 represents the distance between the active aerial target and the origin. Furthermore, 0 is the derivative of 0 . The scalarform of (2) can be written as Where ̇= 0 is assumed considering that fact the main and the pseudo sensors are motionless. Since there are two auxiliary targets, we have a two equations two unknowns system of equations for and . A traditional way to solve this system of equations is to use a grid search, which is a method with high computational complexity, especially when the search area is threedimensional and/or high accuracy for localization is required. Therefore, another method is proposed to solve this problem. In this method, we expand the previous relation by using trigonometric identities, and then obtain the variables explicitly and independently.
As shown in Appendix, we arrive at the equations for and as outlined in the following. Considering the normal flight scenario (no abrupt changes in altitude), we assume to be two-dimensional. We track altitude changes for the aerial platform in subsequent iterations of the algorithm.
To simplify, we use the following substitutions in parts of (4).
After determining the exact value of , we calculate the exact value of the velocity vector angle . 11 At the beginning of this subsection, target coordinates were calculated for each hypothetical altitude, and at the end of this subsection, the velocity vector's magnitude and direction were obtained explicitly and independently. In the next section, we will use a method to solve the problem of target altitude determination.

III. SOLVING FOR TARGET ALTITUDE BY MOTION ANALYSIS
This method, which is based on motion analysis of aerial targets, solves for the target altitude iteratively [19,20]. As we know, Fang's method suffers from considerable levels of error in noisy conditions since it cannot use additional measurements (extra sensors) and cannot solve the problem in overdetermined conditions. Moreover, one cannot localize the target in three-dimensions, using three sensors. Thus, we use motion analysis to increase the localization accuracy of the mentioned method and to resolve the ambiguity of the target altitude.

A. Motion Analysis and Increasing Localization Accuracy
Due to the iterative nature of the algorithm presented in this paper, positioning errors in successive iterations can lead to non-convergence of the algorithm. For this reason and in order to increase the localization accuracy in the first step, we exploit integration of Kalman and UFIR filters to take advantage of accuracy and robustness provided by the two [19]. These operations are performed by the adaptive interacting multiple model (IMM) block which estimates the optimal state vector and covariance matrix [19], as shown in the block diagram of Fig. 3.
Multiple known targets (pseudo-sensors) are fixed in > REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 4 predetermined locations, and one central terrestrial station with receiver channels is used for signal processing to extract TDOA and FDOA values and determine the position and the velocity of the moving aerial target. The prediction of the target state vector at the time step as a function of its previous state is given by is the zero mean process noise, (not necessarily Gaussian) with covariance matrix , and indicates the aerial platform's state which consists of three dimensional position and velocity components of the target as below [ , , , , ,0] Assuming the interval between successive measurements to be short, a constant velocity motion model is considered to estimate the next position of the target. The advantage of using FDOA data in this section is that in the Kalman and UFIR filters' update stage, the magnitude and direction of the velocity vector of the aerial moving target are also corrected. In this case, the measurement vector ( ) becomes , and targets with more complex motion models can be localized using the same assumptions. Therefore, the Kalman and the UFIR filters are applied to the position and speed of the target with a constant speed transfer matrix. The state-transition matrix for this system is. Here the notations × and 0 × indicate × identity matrix and × zero matrix, respectively. Δ refers to the time step and is the white Gaussian noise process.
The observation equation representing from TDOA/FDOA is given as is the white Gaussian noise vector, and is 66 HI      The variances of TDOA and FDOA measurement are denoted by 2 and 2 , respectively.

B. Resolving the Ambiguity of the Target Altitude
As mentioned earlier, we first use the adaptive IMM estimator block to increase the positioning accuracy at any hypothetical altitude, and then, as shown in Fig. 3, the outputs of this block is used along with the prediction of the optimal motion model to calculate for the allocated probability for each altitude.
This method is applied in two levels. After extracting the adaptive IMM estimator block outputs, the prediction for the optimal motion model is made by Time Update block as The function Λ ( ) is defined as the probability density function of the estimated value − ( ), which is assumed to be normal with mean ( ) and covariance ( ). This density function is a normal distribution around the value estimated by the prediction model, as inspired by the method described in [18,19].  In (13), the indices represent the state vectors for each hypothetical altitude (j=1, 2, 3, 4). The main part of the proposed algorithm is the Optimal Fusion and Update of the > REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 5 hypothetical altitude. In this block, the updated probability values assigned to each altitude ( ) are calculated using Λ ( ) and the probabilities assigned to each altitude in the previous iteration. In the first iteration, equal probabilities are considered for the four hypothetical altitudes of the aerial target (for example [0, 5km, 10km, 15km]). In (15), ( ) is the conditional probability of the th state in time step , provided the target was in the th state at the time step − 1. The state here refers to the motion profile at the altitudes assumed for the aerial target. The state vector was defined previously. In fact, as can be seen in the definition of X j , we have not considered the parameters of altitude and speed in the direction of altitude because the altitude is updated in each iteration and is known as the state tag. In the next step, the new inputs of the motion model are calculated according to the conditional probabilities obtained in the previous step.
The parameter ℎ + ( ) in (16) is equal to the updated value of the system state tag (altitude of target) in the previous iteration, which can be the initial value or the output obtained from the previous iteration. As in section , the value ℎ 0 ( − 1) is used to calculate the position and speed of the target.
Finally, according to the values of h(k) calculated in (2) and using the method proposed in section II, the two-dimensional coordinates and the velocity related to the new values of the updated altitudes h(k) are obtained. The pseudocode of the proposed iterative method can be seen below in Algorithm 1. In this pseudocode, first the position of the target and its velocity vector are obtained, and then using the algorithm presented in the previous section, the hypothetical altitudes of the target are updated. This is repeated until the altitude difference between two iterations becomes less than a certain value. This value, which depends on the SNR of the received signal and the target position, is denoted by 0 . The values of ℎ + are updated using   . 16 Eq .

10
End End begin

IV. SIMULATION RESULTS
In this section, numerical simulations are employed to evaluate the performance of the proposed method. For this purpose we assume the following placement for the main and the pseudo sensors. The main sensor is placed at the origin, and the two pseudo sensors are assumed to be placed at (10, 0, 0)km and (10, 10, 0)km respectively. The aerial target is assumed to be situated at (20, 5, ℎ )km, and simulations are performed for different values of ℎ . Also, target speed is set as = = 70 .
We consider the following form for the matrix , so that at each new step, the probability of staying at the same altitude is the same as that of moving to any of the adjacent altitudes for the target. The effect of considering these values for the matrix can be seen in (15).
The initial values of probabilities are also considered as = 1 , in which = ( ) = 4.
First, the target is assumed to be at the altitude of 10 km, and the result is plotted for after repeating all of the above steps once. The result is shown in Fig. 4. As shown in Fig. 4, the target is more likely to be at an altitude of 10km than any other altitude. The algorithm is continued to check the achievable accuracy and improve the resolution of the target altitude report. The result is shown in Fig. 5 after nine repetitions of all steps of the algorithm. As can be seen, the candidate altitude values [0, 5, 10, 15]km in the first iteration (as shown in Fig. 4) change to [9654, 9850, 10045, 10240]m after 9 repetitions in Fig. 5 being more close to the actual altitude of 1okm.
In order to check the results for altitudes other than the default altitudes of [0, 5, 10, 15]km, we assume that the target flies at an altitude of 8.5km, then the result of the algorithm after 9 repetitions is as shown in Fig. 6, which clearly indicates that the probability of being at the altitude of 8560m is more than any other altitude, and this altitude is reported as the altitude of the target by the algorithm. Note that this result is sufficiently close to the actual altitude of 8500m. Fig. 7 shows how the algorithm converges to the altitude of the target at altitudes of 10km and 8.5km. In this Figure and in each iteration, the altitude with the highest probability of is selected (14). As can be seen, in both cases, after some transient states, the altitudes converge from an initial default value to the correct altitude.  Table 1 shows the probability distribution among the four states in each iteration.  In the following, we compare the distance detection error in this study with that of the conventional two-step and linear methods [21]. In these methods, five receivers are used for passive three-dimensional localization of a target using the time and frequency differences of the received signals. In contrast, in the method presented in this study a single receiver along with two pseudo sensors are used. Therefore, to compare the conventional methods with the method presented in this study, only three receivers are considered. Our method is very advantageous to other methods as long as we are limited to three terrestrial receivers. The arrangement of the three receivers to simulate the mentioned methods is such that the first (main) receiver is placed at the origin and the two pseudo sensors are placed in two arbitrary positions. TDOA and FDOA measurements are derived adding Gaussian noise with a mean of zero and diagonal covariance matrices to each parameter's exact value. The diagonal elements of the covariance matrices are considered as 2 = 2 and 2 = 0.1 2 for TDOA and FDOA measurements respectively and the noise terms of TDOA and FDOA measurements are assumed to be uncorrelated. The position and speed of the target are considered as before with an altitude of 8.5km. In the methods that are to be compared with this study, as long as just three receivers are used the aerial target coordinates will be measured in two dimensions only, ignoring the altitude and solving for target distance in two dimensions. Failure to consider the target altitude adds a fixed error to the calculations, which is = √ 2 + 2 + ℎ 2 − √ 2 + 2 = √20 2 + 5 2 + 8.5 2 − √20 2 + 5 2 = 1.68km After 1000 runs of the simulation for each value of noise variance in all three methods, the result of localization accuracy is derived as shown in fig.8. As shown in this figure, the proposed method performs considerably better compared to the two competing methods.

V. CONCLUSION
In this study we considered a single ground based sensor along with two large terrestrial reflectors (pseudo-sensors) to obtain the three dimensional position of an aerial target. To this end, a closed-form solution was presented for the two dimensional position of the target based on TDOA and FDOA equations, and motion analysis was used to derive target altitude despite the fact that all of the main and pseudo stations were placed on the ground.
In addition to being able to derive the three dimensional location of an aerial target based on just one main and few pseudos' ground based sensors, the presented algorithm also has low computational complexity making it attractive for practical applications.
Presented simulations validated performance of the algorithm and showed its superiority to conventional methods in determining position of aerial targets using ground based sensors.

Appendix
This appendix provides a closed and explicit formula for obtaining the magnitude and direction of the velocity vector of the aerial target using (2). To elaborate, first we rewrite equation (2)  Due to the immobility of sensors, the terms ̇ and ̇ are considered zero. As mentioned earlier, represents the input time difference, ̇ represents the input frequency difference, and 0 represents the magnitude of the target position vector, which measures the magnitude ( ) and direction ( ) in (1). The parameter ̇0 is the derivative of the magnitude of the target position vector ( ), which must be written in terms of magnitude and the direction of the velocity vector ( , ). Let's first write the formula for derivative of a vector with respect to time (1/2) .  (21), we achieve an explicit and independent form for each unknown. To do this, the sine and cosine functions of (21)