Distributed Computation for Direct Position ... - Semantic Scholar
1
Distributed Computation for Direct Position Determination Emitter Location Mohammad Pourhomayoun, Student Member, IEEE, Mark L. Fowler, Senior Member, IEEE Abstract—Classical geolocation based on time-difference-ofarrival and frequency-difference-of-arrival (TDOA/FDOA) uses a two stage estimation approach. The single-stage approach Direct Position Determination (DPD) has been proposed to improve accuracy. However, unlike the classical two-stage method, the proposed DPD method does all processing a single node. That is not desirable when computational capabilities are limited and makes the approach non-robust to loss of the central sensor. We develop and assess several DPD variants that address these issues. Index Terms— Source Localization, Emitter Location, Cross Ambiguity Function (CAF), Frequency Difference Of Arrival (FDOA), Time Difference Of Arrival (TDOA), Direct Position Determination (DPD).
I. INTRODUCTION One of the most accurate and common methods for passive radio signal geolocation is based on exploiting the timedifference-of-arrival (TDOA) and frequency-difference-ofarrival (FDOA) between received signals. The classical approach to this method uses two stages to estimate the signal position. In the first stage, TDOA/FDOA are estimated by computing the cross ambiguity function (CAF) and finding the peak of its magnitude surface [1]. In the second stage of the classical method, the TDOA/FDOA estimates are used in statistical processing to locate the emitter [2]. Recently, some new TDOA/FDOA-based methods have been proposed that use only one stage [5], [6], [7], [11], providing improved accuracy. Weiss and Amar [5], [6], [7] introduced a single-stage location method that they call Direct Position Determination (DPD). In independent work, Kay and Vankayalapati [11] developed an identical estimator as part of a generalized likelihood ratio (GLR) detector based on all the received signals. The DPD method centrally processes all the received signals to computes the CAF between all pairings of signals in order to form the so-called Cross Ambiguity Matrices (CAM), whose largest eigenvalues are computed and used to estimate the location. Neither [7] nor [11] addressed the issues of
Manuscript received xxxxxxx. This work was partially supported by the Air Force Research Lab, Rome, NY under contract FA8750-09-2-0068. Approved for Public Release; Distribution Unlimited: 88ABW-2012-0924, 23-Feb-2012 The authors are with Department of Electrical and Computer Engineering, Binghamton Univ., P.O. Box 6000, Binghamton, NY 139026000, Phone: 607-777-6973, [email protected], [email protected].
computation and data transmission for DPD that may set limitations on the practical use DPD; that is the focus here. There are two main issues when DPD is applied in its originally stated centralized formulation. First is the large amount of computations required at a single sensor. This would require at least one sensor to house sufficient computational capability to do the complete DPD processing in a timely fashion and would have a significant impact on the size, weight and power (SWaP) requirement of that sensor. Thus, a decentralized computational approach that strives to spread the computations more uniformly across the sensors is desirable in most practical scenarios. The second issue is the critical mission dependence on the life of the common sensor. In other words, if the common sensor malfunctions or is otherwise incapacitated during the DPD computation stage, no location estimate will be available. Thus, it is desirable to use distributed methods to reduce the critical dependence on any one sensor. In this paper, we propose three different methods aimed at addressing these issues (to varying degrees). The first method applies the Gershgorin theorem at each sensor to allow a distributed computation of an approximation to the DPD location; we call this method “Gershgorin Theorem – DPD” (GT-DPD). The preliminary results of using Gershgorin theorem has been presented in [15]. In this paper, we expand those results by providing a detailed discussion of how data should be shared and used in the sensor network to achieve the best performance with GT-DPD. We also provide a comparison between GT-DPD and the other proposed methods, which is provided in last section of the paper. The second method proposes sharing the CAFs that have been computed throughout the sensor network so that a common sensor can form the CAM and complete the location processing; we call this method “CAF Sharing – DPD” (CSDPD). Furthermore, to reduce the transmission requirement for CAF-sharing, we propose to apply data compression methods tailored to exploit CAF properties, such as the preliminary approaches developed in [8][12][13]. In this paper we expand our earlier results by discussing the communication and computation structures needed to support this method. The third method is a hybrid of the classical method and DPD – called “Hybrid Classical – DPD” (HC-DPD) – that uses some extracted TDOA/FDOA values to aid in the formation of the CAM and then completes location processing using the DPD method of extracting the location from the CAM. This method has not previously been published in our preliminary work. Finally, comparisons will be made to
2 determine the advantages and disadvantages of each method in terms of estimation accuracy, computational load, transmission load and reliability. It should be noted that Fowler and Chen [3] have derived a Fisher information-based method to compress the raw signal data that is shared between sensors for TDOA/FDOA location; such compression could be applied to the original DPD as well as to our proposed distributed methods – however, in our simulation results we do not apply it (to avoid clouding the effects of the methods proposed here).
II. BACKGROUND Suppose that an RF transmitter sends out signal s(t) that is then intercepted by L moving sensors. The equivalent lowpass (ELP) signal (also called the complex envelope) of the received RF signal at the lth sensor is given by
sˆ(t ) is the ELP signal of the transmitted signal s(t), fl is
the Doppler shift,
l
is the delay, and
l
sˆrl (t ) l sˆ(t l )e j 2 fl t wl (t )
wl (t ) is the ELP signal of the noise.
(1)
sˆrm (t ) and
sˆrn (t ) , respectively. The estimate for TDOA/FDOA can be obtained by finding the peak of the magnitude of the CAF; the CAF of the mth and nth signals is given by
CAFmn ( , ) sˆrm (t )sˆ*rn (t )e jt dt .
(2)
Given TDOA/FDOA estimates made from several pairs of signals, classical TDOA/FDOA location methods statistically process them to arrive at a geolocation in an x-y plane (or more generally, in x-y-z space). The remainder of this section provides some needed details about the single-stage localization method DPD and its formulation as developed in [7]. Assume that each sensor collects N time samples sampled with sampling frequency Fs 1/ Ts . Then, as given in [7] the received sampled signal can be written as
[sˆrl (t1 ) , sˆrl (t2 ) , ... , sˆrl (t N )]T
wl
[ wl (t1 ) , wl (t2 ) , ... , wl (t N )]T
is a vector of the N samples of the received signal at
l sensor, sˆ is a vector of the N samples of the transmitted signal, fl is the Doppler shift and Dl is the time sample shift operator by nl ( l / Ts ) samples1. We can write th
Dl Dnl where D is an N N permutation matrix defined as [ D]ij 1 if i j 1 , [ D]0, N 1 1 and [ D]ij 0 otherwise. According to [7] and [11], the estimated transmitter’s position in TDOA/FDOA-based one-stage method is found as follows. Let pi i1 be grid points of possible emitter locations G
Vpi
[ D1HW1H sˆr1 , D2HW2H sˆr 2 , ... ,DL HWL H sˆrL ],
where the delay and Doppler operators correspond to the path between the grid point and the respective sensors 1 to L [7]. For each grid point, then form the L L matrix
Qpi VpHi Vpi
and find its largest eigenvalue, which gives a
set of G largest eigenvalues. The location estimate is the grid point that maximizes the largest eigenvalue, that is
pˆ arg max {max (Q pi )}
Now, suppose that
two sensors m and n receive the ELP signals
sˆrl
sˆrl
e j 2 fltN }
is a complex
amplitude; this is the commonly used narrowband model [10] for RF signal reception in the presence of relative motion between transmitter and receiver. Using a subscript r to indicate the noisy received signal we have
where
diag{e j 2 fl t1 , e j 2 fl t2 , ... ,
in the x-y plane (or more generally in an x-y-z space). For each grid point form the matrix
sˆl (t ) l sˆ(t l )e j 2 fl t where
where
Wl
(3)
pi , i{1, 2, 3, ... , G}
Since all of this processing should be done for each one of the grid points, it is clear that a large amount of computation must be done to find the location and in the centralized formulation of [7] and [11] all of it is done at the common sensor. As shown in [7] and [11], the (k,l)th element of the matrix Q is the value of the CAF between the signals received by sensor k and sensor l and that is why in [11], this matrix is called the Cross Ambiguity Matrix or CAM as stated in the introduction. The elements of the CAM are then given by
[Q pi ]kl [V HV ]kl sˆrk H Wk Dk Dl H Wl H sˆrl [CAFkl ]( i ,i ) where
i and i
(4)
are the corresponding TDOA and FDOA
between sensors k and l with emitter hypothesized as being located at position pi. Note that the diagonal elements in the CAM (k = l) are
sˆrl l Wl Dl sˆ wl
1 As in [7] we assume that the delay can be approximated as being an integer multiple of the sampling instants.
3
[Q pi ]kk [V HV ]kk
L
sˆrk H Wk Dk Dk H Wk H sˆrk sˆrk , 2
(5)
III. GERSHGORIN THEOREM FOR DPD METHOD (GT-DPD) To address the problem that the DPD computation is concentrated at one sensor, in this section we develop a method to distribute the computation and reduce the amount of computation at any one sensor. Since in DPD the emitter location is estimated by computing the largest eigenvalue of the CAM in (4) at each grid point, we can use the eigenvalue approximation result known as the Gershgorin Theorem (GT) [17]. Applying GT for a Hermitian and positive definite matrix to the CAM gives
, max(CAFi ) CAFi
CAF
(6)
max(CAFi )
(7)
ij
j
and then
ˆmax
i
j 1
up the CAFij ’s and then finds the peak of CAFi
which is equal to the energy of the received signals or, equivalently, the peak of the auto ambiguity function (AAF) at the origin. Thus, DPD computes the CAF between all possible pairs of signals (including self-pairing) to populate the CAM. This requires mapping the CAF to the X-Y plane in the same way as used in [4].
i
3- Each sensor i computes CAFi CAFij by adding
can be used as an approximation of the true largest eigenvalue, and that is the essence of our proposed GT-DPD method as follows. Suppose each of the L sensors broadcasts its received signal to all other sensors in the sensor network. Then, each sensor i is able to compute CAFij for all j=1... L and consequently, it is able to compute CAFi j 1 CAFij . Now, L
if we approximate the largest eigenvalue of the CAM by the upper bound on the eigenvalues ( max ˆmax ), then the location estimation will be determined by the point having the largest ˆmax . Here is the outline for the GT-DPD method: 1- Each sensor broadcasts its received signal in the sensor network. 2- Each sensor i computes CAFij ’s in TDOA/FDOA plane and then maps them from TDOA-FDOA plane to X-Y (emitter position) plane (as is done in [4]). The mapping will be done very easily knowing the position and velocity of the sensors and also the grid point position.
(named CAFi , peak ) and its location
( xi , peak , yi , peak )
and then transfers the three numbers
xi , peak , yi , peak and
CAFi , peak to a common sensor (or to all other sensors since there are just three numbers and there is no communication load to transfer them). Note that it is in this step that the Gershgorin Theorem is exploited. 4- According to the (3) and (7), the emitter location estimated is taken as the ( xi , peak , yi , peak ) corresponding to the largest CAFi , peak over all i. Note that in the original DPD method, we need to recompute and form the CAM for each grid point and find the largest eigenvalue of that matrix each time, which requires a very large amount of computation, especially when the number of receiving sensors is large. Moreover, all of these computations would be done at one single sensor. But, the GTDPD method outlined above does not need to form the CAM at all nor does it need to do computationally expensive computations of the largest eigenvalue each time. Thus, in the new method, not only has the costly eigenvalue computation been removed, but also the process is distributed among all receiving sensors. When implementing DPD in a scenario where each sensor has limited computational abilities it is more desirable to distribute the load so that no single sensor is over-burdened than to reduce the total computational complexity across the L sensors. In the original DPD method, the common sensor needs to compute the CAM for each grid point, but since CAM is a Hermitian L L matrix formed by CAFs, it just needs to find all the entries on and above the L ( L 1) main diagonal. This is equivalent to computing 2 CAFs (as non-diagonal elements) and L signal energies (as diagonal elements). Moreover, the common sensor needs to calculate the largest eigenvalue of the matrix for each grid point. On the other hand, in the GT-DPD method, each sensor only needs to compute ( L 1) CAFs and one signal energy. In addition, each sensor has to find the sum and maximum in (7) in order to compute its Gershgorin bound rather than having the centralized sensor compute the largest eigenvalue. Thus, L rather than having one sensor compute ( L 1) CAFs as in 2 the original DPD, in the method proposed here each sensor computes only ( L 1) CAFs; furthermore, each sensor performs a simple Gershgorin estimation rather than a more complicated eigenvalue computation. Although the total computation done in the network is somewhat larger than in DPD, the amount of computation that any one sensor must do is significantly less in the GT-DPD method (the number of CAF computations each sensor does is 2/L of the amount that the commons sensor does in standard DPD) – therefore each sensor will have a lower SWaP requirement.
4 Notice that the distributed nature of the GT-DPD method also improves the survivability of the system since the loss of any one sensor during processing does not totally destroy the ability to provide an estimate since all the other sensors still have their candidate eigenvalue results from using the Gershgorin result. Obviously, losing some of the Gershgorin estimates could reduce the accuracy of estimation, but the rest of the receivers can continue the estimation process with no interruption and still provide a location estimate (although possibly of poorer quality). The GT-DPD development provides insight into two other previously published alternatives to the classical approach [4], [11]. In the GT-DPD method, if we ignore the maximum operator term in (7) and just take ˆmax CAFi for only one arbitrary sensor i, then the results will be equivalent to the socalled CAF-Map method, an ad hoc proposal first presented in [4]. Thus, our Gershgorin-based development provides a sound theoretical rationale for the previously ad hoc development of the CAF-Map method in [4], and also shows that CAF-Map will provide lower accuracy than the full DPD method and the GT-DPD method developed here. Another related approach, called the pair-wise maximum CAF detector, was developed in [11]. It is based on comparing the value of max max CAFij ( , ) with a threshold CAF for only one j
,
arbitrary i as reference sensor. The results in [11] showed that this method also has much lower quality in detection compared to the GLRT detector based on largest eigenvalue; no results were provided in [11] on the location accuracy of the pair-wise maximum CAF method but given its relationship to the GTDPD method it is clear its performance will not be as good. The GT-DPD simulation results for many different cases show that the eigenvalue upper bound is very close to the true largest eigenvalue. However, this approximation lowers the quality of the estimation slightly. Most importantly we see that the GT-DPD method retains DPD’s characteristic of reducing the SNR threshold for accurate location. We examined the effect of the proposed approximation on the estimation accuracy using Monte-Carlo computer simulations (with 500 runs each time). In this simulation, a set of 8 moving sensors and one stationary emitter are placed in a configuration as shown in Figure 1. There exists a cross ambiguity function for all possible pairs of sensors. We used a BPSK signal with carrier frequency 1GHz as the transmitted signal. The sampling frequency is 80 kHz, symbol rate is 10 kb/s and the number of samples is equal to 4096. We run this simulation for various SNRs (-30, -25, -20, -15, -10, -5, 0 dB) with 500 runs for each SNR. Figure 2 shows the effect of the Gershgorinbased eigenvalue approximation on RMS X-Y location error. IV. CAF-SHARING DPD-BASED METHOD (CS-DPD) In this section we develop another approach to implementing Direct Position Determination (DPD) that distributes the CAF computations among all receivers and then applies data compression to the computed CAFs to facilitate transmitting them to a common site where the DPD location processing is completed. It should be noted that in this method
the only approximation that occurs is due to the use of data compression on the CAFs. Thus, in the absence of data compression the accuracy of this method is the same as in the original DPD and for high enough bit rates on the CAF compression the degradation is negligible. However, as the data rates are reduced a drop in performance occurs. In subsection IV-A we first describe the basic method without data compression and then in subsection IV-B specific methods for compressing the computed CAFs are discussed.
Figure 1: Placement of the sensors and the emitter position used for simulation.
(a)
(b)
Figure 2: RMS errors for X and Y versus SNR.
A- Distributed Computation of CAFs As mentioned in section II, the (i,j)th element of the CAM in (3) at grid point p is the value of the CAF between the signals received by sensors i and sensor j at the TDOA/FDOA point corresponding to the grid point p. Thus, sensor i and sensor j can contribute to compute the (i,j)th element of the CAM (as well as the (j,i)th element since the CAM is Hermitian). In other words, each pair of sensors i and j can share their received signals to compute the CAFij in TDOA/FDOA domain. Then the computed CAFs can be transmitted to the common site; subsection IV-B discusses data compression methods for reducing the amount of data needed to be transferred. The common site receives the CAFs and then maps them into the X-Y plane and then uses the results to form the CAM matrix for each location grid point and complete the standard DPD location given in (3).
5 An issue to be considered is how the signal data is shared within the network to allow computation of the needed CAFs. Of course one possibility is for each sensor to transmit its data to all other sensors, but it is not required by the method. Alternatively, as described below, each sensor only needs to receive raw signal data from a subset of sensors. Figure 3 shows a simple case with L = 4 sensors where we see how the sensors share their received signals: the solid arrows show the “one step ahead” sharing while the dashed arrows show “two steps ahead” sharing. The matrix shown represents the CAM (which is Hermitian) and the numbers indicate the indices of the specific CAFs needed. The rectangles next to the sensors show which CAFs are computed at which sensor. The diagonal elements of the CAM are equal to the energy of the received signals at each receiver. Although not shown in Figure 3, the computed CAFs and the individual received signal energies need to be transmitted to a common site (which likely would be one of the sensors but need not be). Figure 4 shows the signal sharing for the case of L = 5 sensors, which also uses one-step ahead and two-steps ahead sharing. Note that for L = 4 sensors only half of the sensors must do twostep-ahead sharing while for L = 5 sensors all five must do two-step-ahead sharing. For both L = 6 and L = 7 we need up to “three steps ahead” sharing and for L = 6 only half of the sensors need to do the highest step sharing while for L = 7 all sensors must do that. In general, if we let S be the maximum “step ahead” needed then for L sensors then we need S L 2 and if L is odd then all sensors must do that largest step ahead while if L is even then only half the sensors do. Thus, in CS-DPD all the sensors contribute (nearly) equally in the process of computing the CAFs.
A14
11 12 13 14 22 23 24 33 34 44
A12
1
2
4
3
A34, A24
A23, A13
Figure 3: Signal sharing for CS-DPD with four sensors. A14, A15
11 12 13 14 15 22 23 24 25 33 34 35 44 45 A45, A35 55
A12, A25
1
2
5
3
A23, A13
The computed CAFs and the sensor received energies are then transmitted to a common site (possibly one of the sensors) for formation of the CAM and the final location processing. Although not required for the method, it is desirable to reduce the amount of data transmission needed for this CAF-sharing step. The next two subsections develop two different data compression methods that exploit specific CAF properties to compress the shared CAFs. B- Exploiting CAF Symmetry and Spectral Properties for Data Compression The CAF is a two-dimensional complex valued function. Thus, we can consider the CAF to be an image (albeit a complex valued image) and apply image compression methods to it. In [8] we used Embedded Zerotree Wavelet (EZW) [9] to independently compress the real part and imaginary part of CAF for use in the CAF-Map method of [4]. Although other wavelet-based data compression methods, like JPEG2000, are also applicable here we focus on EZW because it provides the essential aspects of modern wavelet-based compression without the excessive complexity and flexibility of something like JPEG2000. Compression causes two kinds of distortion on the CAF data. The first is due to the effect of quantization, which can be modeled by an additive noise. The second is the result of the compression throwing away insignificant coefficients of the wavelet transform, which can be roughly modeled as passing the data through a lowpass filter. Note that a typical CAF contains a large main lobe and (usually) various small side lobes. An important aspect for compression of the CAF is that it is a relatively slowly changing function. Much of the fast changing nature comes from the effect of the additive noise in the received signals. Thus, the important part of the CAF to be retained is a spatially low-pass type “image” that should show up in the medium and low frequency parts of the wavelet transform used in EZW. This is particularly good for the EZW method because it will give rise to many zero tree roots in the EZW processing (see [9] for discussion of zero tree roots). Furthermore, because the any discarded high frequency will contain mostly noise, the EZW algorithm will also tend to perform a de-noising operation to the CAF. In addition to these general exploitable “image” characteristics of the CAF, there are some other special properties of the CAF that can be exploited for compression. One of the most important properties of CAF is the symmetry. In [12] we have proved that the noise-free CAF is symmetric around its magnitude peak:
| CAFij ( P , P ) | | CAFij ( P , P ) |
(8)
where ( p , p ) is the peak of the noise-free CAF magnitude. A34, A24
4
Figure 4: Signal sharing for CS-DPD with five sensors
This symmetry can be exploited for data compression. In practice, the noise slightly perturbs this symmetry so we can rewrite (8) as
| CAFij ( P , P ) | | CAFij ( P , P ) | E (9)
6 where E is an error term that is usually small. Thus, using the symmetry property, it is possible to extract the entire CAF magnitude by transmission of only half of the CAF magnitude plus the small residual signal E. We used Monte Carlo simulations (with 500 runs each time) to examine the performance of the proposed CS-DPD method with EZW compression and also with EZW and symmetry. (Note that without compression the CS-DPD method performs exactly like the original DPD method.) In this simulation, a set of 4 moving sensors and one stationary emitter are placed in a configuration as shown in Figure 5. The signals are BPSK with carrier frequency = 1GHz, the sampling frequency = 400 kHz, SNR = -10 dB and the number of signal samples was 65536. The signals were shared as shown in Figure 3 and the six CAFs shown there were computed, compressed and transmitted to a separate common site (i.e., all CAFs were compressed). Two different compression scenarios were examined. The first scenario applied only the EZW algorithm to compress the CAF (labeled “simple compression” in the figures). The second scenario used the EZW algorithm to compress the CAF together with the symmetry property to reduce the amount of transmitted data (labeled “symmetric compression” in the figures). In this case, the EZW algorithm has been applied only on half of the CAF along with the small residual E (as we see in (9)).
In the previous sub-section we treated the CAF as a (complex-valued) image and applied a compression method tailored for images. Alternatively, the CAF can be viewed as a matrix, which motivates using the singular value decomposition (SVD) for compression. For a complexvalued matrix X, the SVD will be r
X i ui viH i 1
where the ui are the left singular vectors, the vi are the right singular vectors, the i are the nonnegative real singular values ordered such that
i i 1
, and r is the number of
non-zero singular values. By truncating the above summation to k r terms, we get a rank-k matrix Xk that approximates X better than any other rank-k matrix in the least-square error sense [13], [14]. Thus, from a data compression viewpoint we can think of Xk as a minimum MSE distortion version of X. The complex valued M N matrix X contains MN complex values or equivalently 2MN real values. However, the truncated SVD representation Xk can be represented with only k(2M +2N+1) real values: we need kM complex values for the
ui ik 1 , kN complex values for the vi ik 1 , and k real values to k represent the singular values i i 1 . Thus, assuming the same number of bits are used to represent the elements in each case, the compression ratio achieved in using Xk rather than X is: MN CR k M N 1 2
Figure 5: Placement of the sensors and the emitter position used for simulation of CS-DPD with EZW-based compression.
The effect of data compression on the RMS error of the location estimation using CS-DPD with compression is shown in Figure 6. In this figure, the RMS error is in meter and the data rate is in bits/element, where an “element” is defined as a real (or imaginary) value of a CAF point. As expected the level of error decreases as the bit rate increases, although the curve is quite flat once above a certain bit rate threshold. Comparing the two curves in each plot shows that the symmetric compression method provides better location accuracy for the same bit rate.
C- Exploiting CAF Rank Properties for SVD-based Data Compression
For example, the compression ratio for a 12832 matrix truncated with k = 1 is 25:1. The CAF usually contains a big main lobe and several small side lobes so that slices at different points we give curves with similar shapes, and this motivates an expectation that a truncation with small k can give an accurate approximation. This idea is supported by results in [16] where it has been shown that time-frequency operators (which include the CAF) have several large singular values followed by a sharp plunge to smaller values and finally decaying asymptotically to zero. Our simulation results also support that most of the information is concentrated in the first few singular vectors/values of the CAF. The effect of noisy signals is to cause small noise-free singular values to be larger even though most of the true CAF information lies concentrated in the first few singular vectors and values. Thus, SVD truncation has the potential to reduce the amount of noise and equivalently increase the effective SNR. To illustrate this, the singular values of a sample 128x32 CAF are illustrated in Figure 7 for two cases: (a) noise-free signals and (b) noisy signals. For the noise-free case, it is clear that there are only 3 to 5 significant singular values, showing that the CAF is very close to a low rank matrix. But for the noisy case the number of significant singular values increases to about 11 or 12. If we truncate the noisy case to have only 3 to 5 terms it is clear that the signal to noise ratio can increase by applying SVD data compression and retaining the first few
7 singular values and discarding the rest.
accurate results than without compression case. This improvement is obtained because of the de-noising property of SVD-based data compression.
(a)
(b)
Figure 7: Singular Values of CAF for two cases: (a) Noise-free signal, (b) Noisy signal
Figure 6: RMS errors for X and Y versus bits/element for CSDPD with EZW-based compression.
We examined the performance of the CS-DPD method using Monte Carlo computer simulations (with 500 runs each time). In this simulation, a set of 4 moving sensors and one stationary emitter are placed in a configuration as shown in Figure 8. In this simulation, the signals are BPSK, the sampling frequency = 20 kHz and the number of signal samples was 4096. As in the EZW case, the signals were shared as shown in Figure 3 and the six CAFs shown there were computed, compressed and transmitted to a separate common site (i.e., all CAFs were compressed). Figure 9 shows the effect of data compression on RMS error of emitter location estimation for the X and Y dimensions. The four curves compare the cases (i) without compression, (ii) SVD-based compression with compression ratio of 25:1, (iii) SVD-based compression with compression ratio of 8:1, and (iv) SVD-based compression with compression ratio of 5:1. The results show that even for a high compression ratio of 25:1 the estimation accuracy is quite close to the case without compression. Surprisingly, at low SNR the case with the compression ratio of 5:1 (and even the case with the compression ratio of 8:1 in some points) yields more
Figure 8: Placement of the sensors and the emitter position used for simulation of CS-DPD with SVD-based compression.
V. HYBRID CLASSICAL AND DPD (HC-DPD) In the classical method of TDOA/FDOA location the CAF computation is typically distributed (and is usually limited to a subset of the total possible CAFs that can be computed) but unlike the CS-DPD method it only transfers the TDOA/FDOA values extracted from the CAFs to the center site rather than complete CAFs as in CS-DPD. Because of the CAF transferal, the DPD method is known to outperform the classical method [7]. In this section, we develop a method taking the advantages of both classical and DPD scenarios and exploits the relationship between different CAFs to reduce the amount of data transmission in the sensor network.
A- Relationships Between CAFs from Various Pairs
8 Suppose that two sensors S1 and S2 receive the ELP signals sˆ1 (t ) and sˆ2 (t ) , respectively. Under the so-called
12 (12 ) CAF12* (12 , 12 ) / A11* (0,0)
(13)
narrowband approximation the ELP model of the received signal for the noise-free case will be [10]:
sˆ1 (t ) 1e jc1 e j1t sˆ(t 1 )
where
i
sˆ2 (t ) 2e jc 2 e j2t sˆ(t 2 )
i is the Doppler shift, i is the and c is the carrier frequency in
is the time delay,
positive-real amplitude
rad/sec [10]. Now, we can write one of the received signals in terms of the other one:
sˆ2 (t )
2 j j j t e e e sˆ1 (t 12 ) 1 c 12
1 12
12
sˆ2 (t ) G12e j12t sˆ1 (t 12 )
(10)
2 j j e e 1
G12
c 12
1 12
where 12 (1 2 ) is the TDOA and 12 (1 2 ) is the FDOA. The CAF between these two signals can be written as being proportional to the Auto Ambiguity function of the first signal as:
CAF12 ( , )
sˆ (t )sˆ (t )e * 2
1
jt
dt
G12* e j12
sˆ (t )sˆ (t * 1
1
12
)e j ( 12 )t dt
which gives
CAF12 ( , ) 12* ( ). A11 ( 12 , 12 )
12 ( ) G12e j
(11)
12
where
A11 ( , )
sˆ (t )sˆ (t )e * 1
1
jt
dt
is
Figure 9: Simulation results showing RMS error for X and Y when we compress CAF using SVD-based compression. The four curves compare the cases (i) without compression, (ii) SVD-based compression with compression ratio of 25:1, (iii) SVD-based compression with co compression ratio of 8:1, and (iv) SVD-based compression with compression ratio of 5:1.
Auto
In general, we can show that:
Ambiguity function (AAF). If there is a third received signal, we can similarly show that,
CAFmn , 1m (1m )e j1m CAF1n 1m , 1m
(14)
CAF23 , sˆ2 (t ) sˆ3* (t )e jt dt
Similarly, we can use 1m (1m ) CAF1*m (1m , 1m ) / A11 (0, 0) to compute the term 1m (1m ) in (14).
G
12
e j12t sˆ1 (t 12 ) sˆ3* (t )e jt dt
12 (12 )e
Further, by putting compute
12
12 (12 ) as,
j12
and
CAF13 12 , 12
12
(12) in (11), we can
B- Using CAF Relationships for Hybrid Classical – DPD Exploiting the (14), it is possible to approximately construct the CAM by transferring only its first row plus the corresponding estimated TDOA/FDOAs (estimated as they are in the classical method). Thus, this distributed method is a
9 hybrid combination of the classical method and the DPD method so we call it the Hybrid Classical – DPD method or HC-DPD. Figure 10 illustrates a simple case of 4 receiving sensors. As we see in Figure 10 (a), sensor 1 (any of the sensors could play this role, we chose sensor 1) distributes its received signal to the other sensors in the network, each of which then computes the CAF between its own signal and the signal from sensor 1. Then, each sensor can transmit the computed CAF to a common site (which can be one of the sensors). In Figure 10 (b), we assumed sensor 2 as the common point. After receiving the CAFs, sensor 2 is able to compute CAF23, having CAF12 (computed by itself) and CAF13 (received from sensor 3) applying (12) and (13). In (12), τ12 and ω12 are estimated by peak of the CAF12 magnitude and * * 12 (12 ) CAF12 (12 , 12 ) / A11 (0, 0) where A11(0,0) is the value of the auto ambiguity function of the received signal by sensor1 at its magnitude peak (that can be also computed by sensor2 or can be transferred to sensor2 as a number). Similarly, sensor 2 is also able to compute CAF24 and CAF34 having CAF12, CAF13 and CAF14 applying (14). Now, sensor 2 has everything needed to approximately form the CAM and do the location estimation.
1
sˆr1
2
1
2
sˆr1
sˆr1
A33 CAF 13
A44 CAF 14
4
3
4
(a)
The one-stage localization method is a major new development in TDOA/FDOA-based emitter location that provides significant improvement in performance at low SNR levels. However, that improvement comes at a cost of significantly more computational complexity. Worse, as DPD was proposed, that complexity is all concentrated at one computing node, which is different from the classical method where the computations are distributed evenly among the sensors. Furthermore, unlike the classical method, the location processing is highly complex (requiring the computation of eigenvalues for each grid point). The large amount of mathematical computations that should be done by a single center point may make some restrictions in processing. Moreover, a complete dependence of the whole process (in both data collection and data processing) on only one single point is also a problem that may reduce the reliability of the system. The method developed here reduce the load of data computation and data transmission by using distributed data computation and processing, applying data compression methods, exploiting the CAF properties, taking advantage of CAF relationships in the sensor network and exploiting some kind of beneficial approximations. The proposed methods are summarized in the following paragraphs as well as in TABLE I: Comparison between DPD and the Proposed Methods, where “+” indicates advantages and “-“ indicates disadvantages. As shown in Table 1, the original DPD method is exceptionally good in accuracy but is exceptionally bad in distributing the computations and is bad in terms of data transmission.
3
(b)
Figure 10: HC-DPD for four receivers.
We examined the proposed method using Monte-Carlo computer simulations (with 500 runs each time). In this simulation, a set of 3 moving sensors and one stationary emitter were placed in a configuration as shown in Figure 11. In this simulation we computed CAF12 and CAF13 directly using Stein’s method [1] and then computed the CAF23 using (12). Then having all of the elements, we formed matrix CAM and apply the DPD method. The signals are BPSK with carrier frequency 1GHz, the sampling frequency is 20 kHz and the number of samples is equal to 4096. Figure 12 shows the RMS error of emitter location estimation for X and Y dimensions compared to original DPD and Classical methods. HC-DPD leads to a large reduction in computation load compared to DPD and large reduction in the amount of data transmission compared to CS-DPD but with the cost of lower accuracy; however, as we can see in Figure 12 it is still much more accurate compared to classical TDOA/FDOA method.
VI. CONCLUSION
Figure 11: Placement of the sensors and the emitter position used for simulation.
The GT-DPD method exploits the simplicity of the Gershgorin theorem to approximately compute the largest eigenvalue without the high cost of exactly computing it. This enables each sensor to locally make its best estimate of the location based on that data it has. These locally-generated estimates are then transmitted to a central location where a final decision is made. Use of the Gershgorin Theorem theorem enables nearly the same accuracy as the original DPD
10 method while being exceptionally good at distributing the computations; data transmission requirements, however, are on par with DPD. In CS-DPD method, we applied the distributed computation idea to divide the computational load among all sensors of the network. In this method, we don’t use any other approximation more than data compression and the quality of this method is the same as the Original DPD for proper bit rates. The extra transmission required is mitigated by data compression so the final transmission requirements are on par with DPD. Finally, in the HC-DPD Method, we used the same idea of Decentralized DPD in addition to exploiting the relationship between different CAFs to eliminate the data redundancy and this idea leads to a large data transmission reduction but at the cost of a reduction in accuracy.
REFERENCES [1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14] [15]
[16]
[17]
Figure 12: RMS errors for X and Y estimates versus SNR
TABLE I: COMPARISON BETWEEN DPD AND THE PROPOSED METHODS
Method
Accuracy
DPD GT-DPD CS-DPD HC-DPD
++ ++ ++ +
Distributed Computation -++ + +
Data Transmission +
Stein, S., “Algorithms for Ambiguity Function Processing,” IEEE Transactions on Acoustics, Speech, and Signal Processing, 29(3), 588599 (1981). D. J. Torrieri, “Statistical theory of passive location system,” IEEE Transactions on Aerospace and Electronic Systems, vol. AES-20, no. 2, March 1984, pp. 183 – 198. Chen, M., Fowler, M. L. “Data Compression for Multiple Parameter Estimation with Application to Emitter Location Systems,” IEEE Transactions on Aerospace and Electronic Systems, vol. 46, no. 1, January 2010, pp. 308 – 322. Hartwell, G., D., “Improved Geo-Spatial Resolution Using a Modified Approach to the Complex Ambiguity Function (CAF)”, Master’s Thesis, Naval Postgraduate School (2005). Weiss, A., “Direct Position Determination of Narrowband Radio Frequency Transmitters,” IEEE Transactions Signal Process. , 11(5), 513-516 (2004). Amar, A., Weiss, A., “Localization of Narrowband Radio Emitters Based on Doppler Frequency Shifts,” IEEE Transactions Signal Process. , 56(11), 5500-5508 (2008). Weiss, A., Amar, A., “Direct Geolocation of Stationary Wide Band Radio Signal Based on Delays and Doppler Shifts,” IEEE Workshop on Statistical Signal Processing, Aug. 31 – Sept. 3, 2009, Cardiff, Wales, UK. M. Pourhomayoun, M. L. Fowler , “Data compression for complex ambiguity function for emitter location”, Proceedings of SPIE Mathematics of Data/Image Coding, Compression, and Encryption with Applications XII, San Diego, Aug 2010. Shapiro, J., M., “Embedded image coding using zerotrees of wavelet coefficients,” IEEE Transactions on Signal Processing, 41(12), 34453462 (1993). R. E. Blahut, “Theory of remote surveillance algorithms, The IMA Volumes in Mathematics and Its Application, Volume 32, Radar and Sonar, 1991. N. Vankayalapati, S. Kay, “Asymptotically Optimal Localization of an Emitter of Low Probability of Intercept Signals Using Distributed Sensors,” to appear in IEEE Transactions on Aerospace and Electronic Systems, vol. 48, no. 1, Jan. 2012. M. Pourhomayoun and M. L. Fowler, “Exploiting Cross Ambiguity Function Properties for Data Compression in Emitter Location Systems,” Conference on Information Sciences and Systems, Johns Hopkins University, March 2011. M. Pourhomayoun and M. L. Fowler, “An SVD Approach for Data Compression in Emitter Location Systems”, Asilomar Conference on Signals, Systems and Computers, Nov 2011. T. K. Moon, W.C. Stirling, Mathematical Methods and Algorithms for Signal Processing, Prentice-Hall, 2000. M. Pourhomayoun and M.L. Fowler, “Sensor Network Distributed Computation for Direct Position Determination,” IEEE Sensor Array and Multichannel Signal Processing Conference (SAM2012), Jun 2012. C. Heil , J. Ramanathan and P. Topiwala, “Asymptotic Singular Value Decay of Time-Frequency Localization Operators”, Wavelet Applications in Signal and Image Processing II, SPIE, 1994. M. Marcus, H. Minc, “A Survey of Matrix Theory and Matrix Inequalities”, Dover Publication, 1992.
Distributed Computation for Direct Position Determination Emitter Location Mohammad Pourhomayoun, Student Member, IEEE, Mark L. Fowler, Senior Member, IEEE Abstract—Classical geolocation based on time-difference-ofarrival and frequency-difference-of-arrival (TDOA/FDOA) uses a two stage estimation approach. The single-stage approach Direct Position Determination (DPD) has been proposed to improve accuracy. However, unlike the classical two-stage method, the proposed DPD method does all processing a single node. That is not desirable when computational capabilities are limited and makes the approach non-robust to loss of the central sensor. We develop and assess several DPD variants that address these issues. Index Terms— Source Localization, Emitter Location, Cross Ambiguity Function (CAF), Frequency Difference Of Arrival (FDOA), Time Difference Of Arrival (TDOA), Direct Position Determination (DPD).
I. INTRODUCTION One of the most accurate and common methods for passive radio signal geolocation is based on exploiting the timedifference-of-arrival (TDOA) and frequency-difference-ofarrival (FDOA) between received signals. The classical approach to this method uses two stages to estimate the signal position. In the first stage, TDOA/FDOA are estimated by computing the cross ambiguity function (CAF) and finding the peak of its magnitude surface [1]. In the second stage of the classical method, the TDOA/FDOA estimates are used in statistical processing to locate the emitter [2]. Recently, some new TDOA/FDOA-based methods have been proposed that use only one stage [5], [6], [7], [11], providing improved accuracy. Weiss and Amar [5], [6], [7] introduced a single-stage location method that they call Direct Position Determination (DPD). In independent work, Kay and Vankayalapati [11] developed an identical estimator as part of a generalized likelihood ratio (GLR) detector based on all the received signals. The DPD method centrally processes all the received signals to computes the CAF between all pairings of signals in order to form the so-called Cross Ambiguity Matrices (CAM), whose largest eigenvalues are computed and used to estimate the location. Neither [7] nor [11] addressed the issues of
Manuscript received xxxxxxx. This work was partially supported by the Air Force Research Lab, Rome, NY under contract FA8750-09-2-0068. Approved for Public Release; Distribution Unlimited: 88ABW-2012-0924, 23-Feb-2012 The authors are with Department of Electrical and Computer Engineering, Binghamton Univ., P.O. Box 6000, Binghamton, NY 139026000, Phone: 607-777-6973, [email protected], [email protected].
computation and data transmission for DPD that may set limitations on the practical use DPD; that is the focus here. There are two main issues when DPD is applied in its originally stated centralized formulation. First is the large amount of computations required at a single sensor. This would require at least one sensor to house sufficient computational capability to do the complete DPD processing in a timely fashion and would have a significant impact on the size, weight and power (SWaP) requirement of that sensor. Thus, a decentralized computational approach that strives to spread the computations more uniformly across the sensors is desirable in most practical scenarios. The second issue is the critical mission dependence on the life of the common sensor. In other words, if the common sensor malfunctions or is otherwise incapacitated during the DPD computation stage, no location estimate will be available. Thus, it is desirable to use distributed methods to reduce the critical dependence on any one sensor. In this paper, we propose three different methods aimed at addressing these issues (to varying degrees). The first method applies the Gershgorin theorem at each sensor to allow a distributed computation of an approximation to the DPD location; we call this method “Gershgorin Theorem – DPD” (GT-DPD). The preliminary results of using Gershgorin theorem has been presented in [15]. In this paper, we expand those results by providing a detailed discussion of how data should be shared and used in the sensor network to achieve the best performance with GT-DPD. We also provide a comparison between GT-DPD and the other proposed methods, which is provided in last section of the paper. The second method proposes sharing the CAFs that have been computed throughout the sensor network so that a common sensor can form the CAM and complete the location processing; we call this method “CAF Sharing – DPD” (CSDPD). Furthermore, to reduce the transmission requirement for CAF-sharing, we propose to apply data compression methods tailored to exploit CAF properties, such as the preliminary approaches developed in [8][12][13]. In this paper we expand our earlier results by discussing the communication and computation structures needed to support this method. The third method is a hybrid of the classical method and DPD – called “Hybrid Classical – DPD” (HC-DPD) – that uses some extracted TDOA/FDOA values to aid in the formation of the CAM and then completes location processing using the DPD method of extracting the location from the CAM. This method has not previously been published in our preliminary work. Finally, comparisons will be made to
2 determine the advantages and disadvantages of each method in terms of estimation accuracy, computational load, transmission load and reliability. It should be noted that Fowler and Chen [3] have derived a Fisher information-based method to compress the raw signal data that is shared between sensors for TDOA/FDOA location; such compression could be applied to the original DPD as well as to our proposed distributed methods – however, in our simulation results we do not apply it (to avoid clouding the effects of the methods proposed here).
II. BACKGROUND Suppose that an RF transmitter sends out signal s(t) that is then intercepted by L moving sensors. The equivalent lowpass (ELP) signal (also called the complex envelope) of the received RF signal at the lth sensor is given by
sˆ(t ) is the ELP signal of the transmitted signal s(t), fl is
the Doppler shift,
l
is the delay, and
l
sˆrl (t ) l sˆ(t l )e j 2 fl t wl (t )
wl (t ) is the ELP signal of the noise.
(1)
sˆrm (t ) and
sˆrn (t ) , respectively. The estimate for TDOA/FDOA can be obtained by finding the peak of the magnitude of the CAF; the CAF of the mth and nth signals is given by
CAFmn ( , ) sˆrm (t )sˆ*rn (t )e jt dt .
(2)
Given TDOA/FDOA estimates made from several pairs of signals, classical TDOA/FDOA location methods statistically process them to arrive at a geolocation in an x-y plane (or more generally, in x-y-z space). The remainder of this section provides some needed details about the single-stage localization method DPD and its formulation as developed in [7]. Assume that each sensor collects N time samples sampled with sampling frequency Fs 1/ Ts . Then, as given in [7] the received sampled signal can be written as
[sˆrl (t1 ) , sˆrl (t2 ) , ... , sˆrl (t N )]T
wl
[ wl (t1 ) , wl (t2 ) , ... , wl (t N )]T
is a vector of the N samples of the received signal at
l sensor, sˆ is a vector of the N samples of the transmitted signal, fl is the Doppler shift and Dl is the time sample shift operator by nl ( l / Ts ) samples1. We can write th
Dl Dnl where D is an N N permutation matrix defined as [ D]ij 1 if i j 1 , [ D]0, N 1 1 and [ D]ij 0 otherwise. According to [7] and [11], the estimated transmitter’s position in TDOA/FDOA-based one-stage method is found as follows. Let pi i1 be grid points of possible emitter locations G
Vpi
[ D1HW1H sˆr1 , D2HW2H sˆr 2 , ... ,DL HWL H sˆrL ],
where the delay and Doppler operators correspond to the path between the grid point and the respective sensors 1 to L [7]. For each grid point, then form the L L matrix
Qpi VpHi Vpi
and find its largest eigenvalue, which gives a
set of G largest eigenvalues. The location estimate is the grid point that maximizes the largest eigenvalue, that is
pˆ arg max {max (Q pi )}
Now, suppose that
two sensors m and n receive the ELP signals
sˆrl
sˆrl
e j 2 fltN }
is a complex
amplitude; this is the commonly used narrowband model [10] for RF signal reception in the presence of relative motion between transmitter and receiver. Using a subscript r to indicate the noisy received signal we have
where
diag{e j 2 fl t1 , e j 2 fl t2 , ... ,
in the x-y plane (or more generally in an x-y-z space). For each grid point form the matrix
sˆl (t ) l sˆ(t l )e j 2 fl t where
where
Wl
(3)
pi , i{1, 2, 3, ... , G}
Since all of this processing should be done for each one of the grid points, it is clear that a large amount of computation must be done to find the location and in the centralized formulation of [7] and [11] all of it is done at the common sensor. As shown in [7] and [11], the (k,l)th element of the matrix Q is the value of the CAF between the signals received by sensor k and sensor l and that is why in [11], this matrix is called the Cross Ambiguity Matrix or CAM as stated in the introduction. The elements of the CAM are then given by
[Q pi ]kl [V HV ]kl sˆrk H Wk Dk Dl H Wl H sˆrl [CAFkl ]( i ,i ) where
i and i
(4)
are the corresponding TDOA and FDOA
between sensors k and l with emitter hypothesized as being located at position pi. Note that the diagonal elements in the CAM (k = l) are
sˆrl l Wl Dl sˆ wl
1 As in [7] we assume that the delay can be approximated as being an integer multiple of the sampling instants.
3
[Q pi ]kk [V HV ]kk
L
sˆrk H Wk Dk Dk H Wk H sˆrk sˆrk , 2
(5)
III. GERSHGORIN THEOREM FOR DPD METHOD (GT-DPD) To address the problem that the DPD computation is concentrated at one sensor, in this section we develop a method to distribute the computation and reduce the amount of computation at any one sensor. Since in DPD the emitter location is estimated by computing the largest eigenvalue of the CAM in (4) at each grid point, we can use the eigenvalue approximation result known as the Gershgorin Theorem (GT) [17]. Applying GT for a Hermitian and positive definite matrix to the CAM gives
, max(CAFi ) CAFi
CAF
(6)
max(CAFi )
(7)
ij
j
and then
ˆmax
i
j 1
up the CAFij ’s and then finds the peak of CAFi
which is equal to the energy of the received signals or, equivalently, the peak of the auto ambiguity function (AAF) at the origin. Thus, DPD computes the CAF between all possible pairs of signals (including self-pairing) to populate the CAM. This requires mapping the CAF to the X-Y plane in the same way as used in [4].
i
3- Each sensor i computes CAFi CAFij by adding
can be used as an approximation of the true largest eigenvalue, and that is the essence of our proposed GT-DPD method as follows. Suppose each of the L sensors broadcasts its received signal to all other sensors in the sensor network. Then, each sensor i is able to compute CAFij for all j=1... L and consequently, it is able to compute CAFi j 1 CAFij . Now, L
if we approximate the largest eigenvalue of the CAM by the upper bound on the eigenvalues ( max ˆmax ), then the location estimation will be determined by the point having the largest ˆmax . Here is the outline for the GT-DPD method: 1- Each sensor broadcasts its received signal in the sensor network. 2- Each sensor i computes CAFij ’s in TDOA/FDOA plane and then maps them from TDOA-FDOA plane to X-Y (emitter position) plane (as is done in [4]). The mapping will be done very easily knowing the position and velocity of the sensors and also the grid point position.
(named CAFi , peak ) and its location
( xi , peak , yi , peak )
and then transfers the three numbers
xi , peak , yi , peak and
CAFi , peak to a common sensor (or to all other sensors since there are just three numbers and there is no communication load to transfer them). Note that it is in this step that the Gershgorin Theorem is exploited. 4- According to the (3) and (7), the emitter location estimated is taken as the ( xi , peak , yi , peak ) corresponding to the largest CAFi , peak over all i. Note that in the original DPD method, we need to recompute and form the CAM for each grid point and find the largest eigenvalue of that matrix each time, which requires a very large amount of computation, especially when the number of receiving sensors is large. Moreover, all of these computations would be done at one single sensor. But, the GTDPD method outlined above does not need to form the CAM at all nor does it need to do computationally expensive computations of the largest eigenvalue each time. Thus, in the new method, not only has the costly eigenvalue computation been removed, but also the process is distributed among all receiving sensors. When implementing DPD in a scenario where each sensor has limited computational abilities it is more desirable to distribute the load so that no single sensor is over-burdened than to reduce the total computational complexity across the L sensors. In the original DPD method, the common sensor needs to compute the CAM for each grid point, but since CAM is a Hermitian L L matrix formed by CAFs, it just needs to find all the entries on and above the L ( L 1) main diagonal. This is equivalent to computing 2 CAFs (as non-diagonal elements) and L signal energies (as diagonal elements). Moreover, the common sensor needs to calculate the largest eigenvalue of the matrix for each grid point. On the other hand, in the GT-DPD method, each sensor only needs to compute ( L 1) CAFs and one signal energy. In addition, each sensor has to find the sum and maximum in (7) in order to compute its Gershgorin bound rather than having the centralized sensor compute the largest eigenvalue. Thus, L rather than having one sensor compute ( L 1) CAFs as in 2 the original DPD, in the method proposed here each sensor computes only ( L 1) CAFs; furthermore, each sensor performs a simple Gershgorin estimation rather than a more complicated eigenvalue computation. Although the total computation done in the network is somewhat larger than in DPD, the amount of computation that any one sensor must do is significantly less in the GT-DPD method (the number of CAF computations each sensor does is 2/L of the amount that the commons sensor does in standard DPD) – therefore each sensor will have a lower SWaP requirement.
4 Notice that the distributed nature of the GT-DPD method also improves the survivability of the system since the loss of any one sensor during processing does not totally destroy the ability to provide an estimate since all the other sensors still have their candidate eigenvalue results from using the Gershgorin result. Obviously, losing some of the Gershgorin estimates could reduce the accuracy of estimation, but the rest of the receivers can continue the estimation process with no interruption and still provide a location estimate (although possibly of poorer quality). The GT-DPD development provides insight into two other previously published alternatives to the classical approach [4], [11]. In the GT-DPD method, if we ignore the maximum operator term in (7) and just take ˆmax CAFi for only one arbitrary sensor i, then the results will be equivalent to the socalled CAF-Map method, an ad hoc proposal first presented in [4]. Thus, our Gershgorin-based development provides a sound theoretical rationale for the previously ad hoc development of the CAF-Map method in [4], and also shows that CAF-Map will provide lower accuracy than the full DPD method and the GT-DPD method developed here. Another related approach, called the pair-wise maximum CAF detector, was developed in [11]. It is based on comparing the value of max max CAFij ( , ) with a threshold CAF for only one j
,
arbitrary i as reference sensor. The results in [11] showed that this method also has much lower quality in detection compared to the GLRT detector based on largest eigenvalue; no results were provided in [11] on the location accuracy of the pair-wise maximum CAF method but given its relationship to the GTDPD method it is clear its performance will not be as good. The GT-DPD simulation results for many different cases show that the eigenvalue upper bound is very close to the true largest eigenvalue. However, this approximation lowers the quality of the estimation slightly. Most importantly we see that the GT-DPD method retains DPD’s characteristic of reducing the SNR threshold for accurate location. We examined the effect of the proposed approximation on the estimation accuracy using Monte-Carlo computer simulations (with 500 runs each time). In this simulation, a set of 8 moving sensors and one stationary emitter are placed in a configuration as shown in Figure 1. There exists a cross ambiguity function for all possible pairs of sensors. We used a BPSK signal with carrier frequency 1GHz as the transmitted signal. The sampling frequency is 80 kHz, symbol rate is 10 kb/s and the number of samples is equal to 4096. We run this simulation for various SNRs (-30, -25, -20, -15, -10, -5, 0 dB) with 500 runs for each SNR. Figure 2 shows the effect of the Gershgorinbased eigenvalue approximation on RMS X-Y location error. IV. CAF-SHARING DPD-BASED METHOD (CS-DPD) In this section we develop another approach to implementing Direct Position Determination (DPD) that distributes the CAF computations among all receivers and then applies data compression to the computed CAFs to facilitate transmitting them to a common site where the DPD location processing is completed. It should be noted that in this method
the only approximation that occurs is due to the use of data compression on the CAFs. Thus, in the absence of data compression the accuracy of this method is the same as in the original DPD and for high enough bit rates on the CAF compression the degradation is negligible. However, as the data rates are reduced a drop in performance occurs. In subsection IV-A we first describe the basic method without data compression and then in subsection IV-B specific methods for compressing the computed CAFs are discussed.
Figure 1: Placement of the sensors and the emitter position used for simulation.
(a)
(b)
Figure 2: RMS errors for X and Y versus SNR.
A- Distributed Computation of CAFs As mentioned in section II, the (i,j)th element of the CAM in (3) at grid point p is the value of the CAF between the signals received by sensors i and sensor j at the TDOA/FDOA point corresponding to the grid point p. Thus, sensor i and sensor j can contribute to compute the (i,j)th element of the CAM (as well as the (j,i)th element since the CAM is Hermitian). In other words, each pair of sensors i and j can share their received signals to compute the CAFij in TDOA/FDOA domain. Then the computed CAFs can be transmitted to the common site; subsection IV-B discusses data compression methods for reducing the amount of data needed to be transferred. The common site receives the CAFs and then maps them into the X-Y plane and then uses the results to form the CAM matrix for each location grid point and complete the standard DPD location given in (3).
5 An issue to be considered is how the signal data is shared within the network to allow computation of the needed CAFs. Of course one possibility is for each sensor to transmit its data to all other sensors, but it is not required by the method. Alternatively, as described below, each sensor only needs to receive raw signal data from a subset of sensors. Figure 3 shows a simple case with L = 4 sensors where we see how the sensors share their received signals: the solid arrows show the “one step ahead” sharing while the dashed arrows show “two steps ahead” sharing. The matrix shown represents the CAM (which is Hermitian) and the numbers indicate the indices of the specific CAFs needed. The rectangles next to the sensors show which CAFs are computed at which sensor. The diagonal elements of the CAM are equal to the energy of the received signals at each receiver. Although not shown in Figure 3, the computed CAFs and the individual received signal energies need to be transmitted to a common site (which likely would be one of the sensors but need not be). Figure 4 shows the signal sharing for the case of L = 5 sensors, which also uses one-step ahead and two-steps ahead sharing. Note that for L = 4 sensors only half of the sensors must do twostep-ahead sharing while for L = 5 sensors all five must do two-step-ahead sharing. For both L = 6 and L = 7 we need up to “three steps ahead” sharing and for L = 6 only half of the sensors need to do the highest step sharing while for L = 7 all sensors must do that. In general, if we let S be the maximum “step ahead” needed then for L sensors then we need S L 2 and if L is odd then all sensors must do that largest step ahead while if L is even then only half the sensors do. Thus, in CS-DPD all the sensors contribute (nearly) equally in the process of computing the CAFs.
A14
11 12 13 14 22 23 24 33 34 44
A12
1
2
4
3
A34, A24
A23, A13
Figure 3: Signal sharing for CS-DPD with four sensors. A14, A15
11 12 13 14 15 22 23 24 25 33 34 35 44 45 A45, A35 55
A12, A25
1
2
5
3
A23, A13
The computed CAFs and the sensor received energies are then transmitted to a common site (possibly one of the sensors) for formation of the CAM and the final location processing. Although not required for the method, it is desirable to reduce the amount of data transmission needed for this CAF-sharing step. The next two subsections develop two different data compression methods that exploit specific CAF properties to compress the shared CAFs. B- Exploiting CAF Symmetry and Spectral Properties for Data Compression The CAF is a two-dimensional complex valued function. Thus, we can consider the CAF to be an image (albeit a complex valued image) and apply image compression methods to it. In [8] we used Embedded Zerotree Wavelet (EZW) [9] to independently compress the real part and imaginary part of CAF for use in the CAF-Map method of [4]. Although other wavelet-based data compression methods, like JPEG2000, are also applicable here we focus on EZW because it provides the essential aspects of modern wavelet-based compression without the excessive complexity and flexibility of something like JPEG2000. Compression causes two kinds of distortion on the CAF data. The first is due to the effect of quantization, which can be modeled by an additive noise. The second is the result of the compression throwing away insignificant coefficients of the wavelet transform, which can be roughly modeled as passing the data through a lowpass filter. Note that a typical CAF contains a large main lobe and (usually) various small side lobes. An important aspect for compression of the CAF is that it is a relatively slowly changing function. Much of the fast changing nature comes from the effect of the additive noise in the received signals. Thus, the important part of the CAF to be retained is a spatially low-pass type “image” that should show up in the medium and low frequency parts of the wavelet transform used in EZW. This is particularly good for the EZW method because it will give rise to many zero tree roots in the EZW processing (see [9] for discussion of zero tree roots). Furthermore, because the any discarded high frequency will contain mostly noise, the EZW algorithm will also tend to perform a de-noising operation to the CAF. In addition to these general exploitable “image” characteristics of the CAF, there are some other special properties of the CAF that can be exploited for compression. One of the most important properties of CAF is the symmetry. In [12] we have proved that the noise-free CAF is symmetric around its magnitude peak:
| CAFij ( P , P ) | | CAFij ( P , P ) |
(8)
where ( p , p ) is the peak of the noise-free CAF magnitude. A34, A24
4
Figure 4: Signal sharing for CS-DPD with five sensors
This symmetry can be exploited for data compression. In practice, the noise slightly perturbs this symmetry so we can rewrite (8) as
| CAFij ( P , P ) | | CAFij ( P , P ) | E (9)
6 where E is an error term that is usually small. Thus, using the symmetry property, it is possible to extract the entire CAF magnitude by transmission of only half of the CAF magnitude plus the small residual signal E. We used Monte Carlo simulations (with 500 runs each time) to examine the performance of the proposed CS-DPD method with EZW compression and also with EZW and symmetry. (Note that without compression the CS-DPD method performs exactly like the original DPD method.) In this simulation, a set of 4 moving sensors and one stationary emitter are placed in a configuration as shown in Figure 5. The signals are BPSK with carrier frequency = 1GHz, the sampling frequency = 400 kHz, SNR = -10 dB and the number of signal samples was 65536. The signals were shared as shown in Figure 3 and the six CAFs shown there were computed, compressed and transmitted to a separate common site (i.e., all CAFs were compressed). Two different compression scenarios were examined. The first scenario applied only the EZW algorithm to compress the CAF (labeled “simple compression” in the figures). The second scenario used the EZW algorithm to compress the CAF together with the symmetry property to reduce the amount of transmitted data (labeled “symmetric compression” in the figures). In this case, the EZW algorithm has been applied only on half of the CAF along with the small residual E (as we see in (9)).
In the previous sub-section we treated the CAF as a (complex-valued) image and applied a compression method tailored for images. Alternatively, the CAF can be viewed as a matrix, which motivates using the singular value decomposition (SVD) for compression. For a complexvalued matrix X, the SVD will be r
X i ui viH i 1
where the ui are the left singular vectors, the vi are the right singular vectors, the i are the nonnegative real singular values ordered such that
i i 1
, and r is the number of
non-zero singular values. By truncating the above summation to k r terms, we get a rank-k matrix Xk that approximates X better than any other rank-k matrix in the least-square error sense [13], [14]. Thus, from a data compression viewpoint we can think of Xk as a minimum MSE distortion version of X. The complex valued M N matrix X contains MN complex values or equivalently 2MN real values. However, the truncated SVD representation Xk can be represented with only k(2M +2N+1) real values: we need kM complex values for the
ui ik 1 , kN complex values for the vi ik 1 , and k real values to k represent the singular values i i 1 . Thus, assuming the same number of bits are used to represent the elements in each case, the compression ratio achieved in using Xk rather than X is: MN CR k M N 1 2
Figure 5: Placement of the sensors and the emitter position used for simulation of CS-DPD with EZW-based compression.
The effect of data compression on the RMS error of the location estimation using CS-DPD with compression is shown in Figure 6. In this figure, the RMS error is in meter and the data rate is in bits/element, where an “element” is defined as a real (or imaginary) value of a CAF point. As expected the level of error decreases as the bit rate increases, although the curve is quite flat once above a certain bit rate threshold. Comparing the two curves in each plot shows that the symmetric compression method provides better location accuracy for the same bit rate.
C- Exploiting CAF Rank Properties for SVD-based Data Compression
For example, the compression ratio for a 12832 matrix truncated with k = 1 is 25:1. The CAF usually contains a big main lobe and several small side lobes so that slices at different points we give curves with similar shapes, and this motivates an expectation that a truncation with small k can give an accurate approximation. This idea is supported by results in [16] where it has been shown that time-frequency operators (which include the CAF) have several large singular values followed by a sharp plunge to smaller values and finally decaying asymptotically to zero. Our simulation results also support that most of the information is concentrated in the first few singular vectors/values of the CAF. The effect of noisy signals is to cause small noise-free singular values to be larger even though most of the true CAF information lies concentrated in the first few singular vectors and values. Thus, SVD truncation has the potential to reduce the amount of noise and equivalently increase the effective SNR. To illustrate this, the singular values of a sample 128x32 CAF are illustrated in Figure 7 for two cases: (a) noise-free signals and (b) noisy signals. For the noise-free case, it is clear that there are only 3 to 5 significant singular values, showing that the CAF is very close to a low rank matrix. But for the noisy case the number of significant singular values increases to about 11 or 12. If we truncate the noisy case to have only 3 to 5 terms it is clear that the signal to noise ratio can increase by applying SVD data compression and retaining the first few
7 singular values and discarding the rest.
accurate results than without compression case. This improvement is obtained because of the de-noising property of SVD-based data compression.
(a)
(b)
Figure 7: Singular Values of CAF for two cases: (a) Noise-free signal, (b) Noisy signal
Figure 6: RMS errors for X and Y versus bits/element for CSDPD with EZW-based compression.
We examined the performance of the CS-DPD method using Monte Carlo computer simulations (with 500 runs each time). In this simulation, a set of 4 moving sensors and one stationary emitter are placed in a configuration as shown in Figure 8. In this simulation, the signals are BPSK, the sampling frequency = 20 kHz and the number of signal samples was 4096. As in the EZW case, the signals were shared as shown in Figure 3 and the six CAFs shown there were computed, compressed and transmitted to a separate common site (i.e., all CAFs were compressed). Figure 9 shows the effect of data compression on RMS error of emitter location estimation for the X and Y dimensions. The four curves compare the cases (i) without compression, (ii) SVD-based compression with compression ratio of 25:1, (iii) SVD-based compression with compression ratio of 8:1, and (iv) SVD-based compression with compression ratio of 5:1. The results show that even for a high compression ratio of 25:1 the estimation accuracy is quite close to the case without compression. Surprisingly, at low SNR the case with the compression ratio of 5:1 (and even the case with the compression ratio of 8:1 in some points) yields more
Figure 8: Placement of the sensors and the emitter position used for simulation of CS-DPD with SVD-based compression.
V. HYBRID CLASSICAL AND DPD (HC-DPD) In the classical method of TDOA/FDOA location the CAF computation is typically distributed (and is usually limited to a subset of the total possible CAFs that can be computed) but unlike the CS-DPD method it only transfers the TDOA/FDOA values extracted from the CAFs to the center site rather than complete CAFs as in CS-DPD. Because of the CAF transferal, the DPD method is known to outperform the classical method [7]. In this section, we develop a method taking the advantages of both classical and DPD scenarios and exploits the relationship between different CAFs to reduce the amount of data transmission in the sensor network.
A- Relationships Between CAFs from Various Pairs
8 Suppose that two sensors S1 and S2 receive the ELP signals sˆ1 (t ) and sˆ2 (t ) , respectively. Under the so-called
12 (12 ) CAF12* (12 , 12 ) / A11* (0,0)
(13)
narrowband approximation the ELP model of the received signal for the noise-free case will be [10]:
sˆ1 (t ) 1e jc1 e j1t sˆ(t 1 )
where
i
sˆ2 (t ) 2e jc 2 e j2t sˆ(t 2 )
i is the Doppler shift, i is the and c is the carrier frequency in
is the time delay,
positive-real amplitude
rad/sec [10]. Now, we can write one of the received signals in terms of the other one:
sˆ2 (t )
2 j j j t e e e sˆ1 (t 12 ) 1 c 12
1 12
12
sˆ2 (t ) G12e j12t sˆ1 (t 12 )
(10)
2 j j e e 1
G12
c 12
1 12
where 12 (1 2 ) is the TDOA and 12 (1 2 ) is the FDOA. The CAF between these two signals can be written as being proportional to the Auto Ambiguity function of the first signal as:
CAF12 ( , )
sˆ (t )sˆ (t )e * 2
1
jt
dt
G12* e j12
sˆ (t )sˆ (t * 1
1
12
)e j ( 12 )t dt
which gives
CAF12 ( , ) 12* ( ). A11 ( 12 , 12 )
12 ( ) G12e j
(11)
12
where
A11 ( , )
sˆ (t )sˆ (t )e * 1
1
jt
dt
is
Figure 9: Simulation results showing RMS error for X and Y when we compress CAF using SVD-based compression. The four curves compare the cases (i) without compression, (ii) SVD-based compression with compression ratio of 25:1, (iii) SVD-based compression with co compression ratio of 8:1, and (iv) SVD-based compression with compression ratio of 5:1.
Auto
In general, we can show that:
Ambiguity function (AAF). If there is a third received signal, we can similarly show that,
CAFmn , 1m (1m )e j1m CAF1n 1m , 1m
(14)
CAF23 , sˆ2 (t ) sˆ3* (t )e jt dt
Similarly, we can use 1m (1m ) CAF1*m (1m , 1m ) / A11 (0, 0) to compute the term 1m (1m ) in (14).
G
12
e j12t sˆ1 (t 12 ) sˆ3* (t )e jt dt
12 (12 )e
Further, by putting compute
12
12 (12 ) as,
j12
and
CAF13 12 , 12
12
(12) in (11), we can
B- Using CAF Relationships for Hybrid Classical – DPD Exploiting the (14), it is possible to approximately construct the CAM by transferring only its first row plus the corresponding estimated TDOA/FDOAs (estimated as they are in the classical method). Thus, this distributed method is a
9 hybrid combination of the classical method and the DPD method so we call it the Hybrid Classical – DPD method or HC-DPD. Figure 10 illustrates a simple case of 4 receiving sensors. As we see in Figure 10 (a), sensor 1 (any of the sensors could play this role, we chose sensor 1) distributes its received signal to the other sensors in the network, each of which then computes the CAF between its own signal and the signal from sensor 1. Then, each sensor can transmit the computed CAF to a common site (which can be one of the sensors). In Figure 10 (b), we assumed sensor 2 as the common point. After receiving the CAFs, sensor 2 is able to compute CAF23, having CAF12 (computed by itself) and CAF13 (received from sensor 3) applying (12) and (13). In (12), τ12 and ω12 are estimated by peak of the CAF12 magnitude and * * 12 (12 ) CAF12 (12 , 12 ) / A11 (0, 0) where A11(0,0) is the value of the auto ambiguity function of the received signal by sensor1 at its magnitude peak (that can be also computed by sensor2 or can be transferred to sensor2 as a number). Similarly, sensor 2 is also able to compute CAF24 and CAF34 having CAF12, CAF13 and CAF14 applying (14). Now, sensor 2 has everything needed to approximately form the CAM and do the location estimation.
1
sˆr1
2
1
2
sˆr1
sˆr1
A33 CAF 13
A44 CAF 14
4
3
4
(a)
The one-stage localization method is a major new development in TDOA/FDOA-based emitter location that provides significant improvement in performance at low SNR levels. However, that improvement comes at a cost of significantly more computational complexity. Worse, as DPD was proposed, that complexity is all concentrated at one computing node, which is different from the classical method where the computations are distributed evenly among the sensors. Furthermore, unlike the classical method, the location processing is highly complex (requiring the computation of eigenvalues for each grid point). The large amount of mathematical computations that should be done by a single center point may make some restrictions in processing. Moreover, a complete dependence of the whole process (in both data collection and data processing) on only one single point is also a problem that may reduce the reliability of the system. The method developed here reduce the load of data computation and data transmission by using distributed data computation and processing, applying data compression methods, exploiting the CAF properties, taking advantage of CAF relationships in the sensor network and exploiting some kind of beneficial approximations. The proposed methods are summarized in the following paragraphs as well as in TABLE I: Comparison between DPD and the Proposed Methods, where “+” indicates advantages and “-“ indicates disadvantages. As shown in Table 1, the original DPD method is exceptionally good in accuracy but is exceptionally bad in distributing the computations and is bad in terms of data transmission.
3
(b)
Figure 10: HC-DPD for four receivers.
We examined the proposed method using Monte-Carlo computer simulations (with 500 runs each time). In this simulation, a set of 3 moving sensors and one stationary emitter were placed in a configuration as shown in Figure 11. In this simulation we computed CAF12 and CAF13 directly using Stein’s method [1] and then computed the CAF23 using (12). Then having all of the elements, we formed matrix CAM and apply the DPD method. The signals are BPSK with carrier frequency 1GHz, the sampling frequency is 20 kHz and the number of samples is equal to 4096. Figure 12 shows the RMS error of emitter location estimation for X and Y dimensions compared to original DPD and Classical methods. HC-DPD leads to a large reduction in computation load compared to DPD and large reduction in the amount of data transmission compared to CS-DPD but with the cost of lower accuracy; however, as we can see in Figure 12 it is still much more accurate compared to classical TDOA/FDOA method.
VI. CONCLUSION
Figure 11: Placement of the sensors and the emitter position used for simulation.
The GT-DPD method exploits the simplicity of the Gershgorin theorem to approximately compute the largest eigenvalue without the high cost of exactly computing it. This enables each sensor to locally make its best estimate of the location based on that data it has. These locally-generated estimates are then transmitted to a central location where a final decision is made. Use of the Gershgorin Theorem theorem enables nearly the same accuracy as the original DPD
10 method while being exceptionally good at distributing the computations; data transmission requirements, however, are on par with DPD. In CS-DPD method, we applied the distributed computation idea to divide the computational load among all sensors of the network. In this method, we don’t use any other approximation more than data compression and the quality of this method is the same as the Original DPD for proper bit rates. The extra transmission required is mitigated by data compression so the final transmission requirements are on par with DPD. Finally, in the HC-DPD Method, we used the same idea of Decentralized DPD in addition to exploiting the relationship between different CAFs to eliminate the data redundancy and this idea leads to a large data transmission reduction but at the cost of a reduction in accuracy.
REFERENCES [1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14] [15]
[16]
[17]
Figure 12: RMS errors for X and Y estimates versus SNR
TABLE I: COMPARISON BETWEEN DPD AND THE PROPOSED METHODS
Method
Accuracy
DPD GT-DPD CS-DPD HC-DPD
++ ++ ++ +
Distributed Computation -++ + +
Data Transmission +
Stein, S., “Algorithms for Ambiguity Function Processing,” IEEE Transactions on Acoustics, Speech, and Signal Processing, 29(3), 588599 (1981). D. J. Torrieri, “Statistical theory of passive location system,” IEEE Transactions on Aerospace and Electronic Systems, vol. AES-20, no. 2, March 1984, pp. 183 – 198. Chen, M., Fowler, M. L. “Data Compression for Multiple Parameter Estimation with Application to Emitter Location Systems,” IEEE Transactions on Aerospace and Electronic Systems, vol. 46, no. 1, January 2010, pp. 308 – 322. Hartwell, G., D., “Improved Geo-Spatial Resolution Using a Modified Approach to the Complex Ambiguity Function (CAF)”, Master’s Thesis, Naval Postgraduate School (2005). Weiss, A., “Direct Position Determination of Narrowband Radio Frequency Transmitters,” IEEE Transactions Signal Process. , 11(5), 513-516 (2004). Amar, A., Weiss, A., “Localization of Narrowband Radio Emitters Based on Doppler Frequency Shifts,” IEEE Transactions Signal Process. , 56(11), 5500-5508 (2008). Weiss, A., Amar, A., “Direct Geolocation of Stationary Wide Band Radio Signal Based on Delays and Doppler Shifts,” IEEE Workshop on Statistical Signal Processing, Aug. 31 – Sept. 3, 2009, Cardiff, Wales, UK. M. Pourhomayoun, M. L. Fowler , “Data compression for complex ambiguity function for emitter location”, Proceedings of SPIE Mathematics of Data/Image Coding, Compression, and Encryption with Applications XII, San Diego, Aug 2010. Shapiro, J., M., “Embedded image coding using zerotrees of wavelet coefficients,” IEEE Transactions on Signal Processing, 41(12), 34453462 (1993). R. E. Blahut, “Theory of remote surveillance algorithms, The IMA Volumes in Mathematics and Its Application, Volume 32, Radar and Sonar, 1991. N. Vankayalapati, S. Kay, “Asymptotically Optimal Localization of an Emitter of Low Probability of Intercept Signals Using Distributed Sensors,” to appear in IEEE Transactions on Aerospace and Electronic Systems, vol. 48, no. 1, Jan. 2012. M. Pourhomayoun and M. L. Fowler, “Exploiting Cross Ambiguity Function Properties for Data Compression in Emitter Location Systems,” Conference on Information Sciences and Systems, Johns Hopkins University, March 2011. M. Pourhomayoun and M. L. Fowler, “An SVD Approach for Data Compression in Emitter Location Systems”, Asilomar Conference on Signals, Systems and Computers, Nov 2011. T. K. Moon, W.C. Stirling, Mathematical Methods and Algorithms for Signal Processing, Prentice-Hall, 2000. M. Pourhomayoun and M.L. Fowler, “Sensor Network Distributed Computation for Direct Position Determination,” IEEE Sensor Array and Multichannel Signal Processing Conference (SAM2012), Jun 2012. C. Heil , J. Ramanathan and P. Topiwala, “Asymptotic Singular Value Decay of Time-Frequency Localization Operators”, Wavelet Applications in Signal and Image Processing II, SPIE, 1994. M. Marcus, H. Minc, “A Survey of Matrix Theory and Matrix Inequalities”, Dover Publication, 1992.
