A Novel Virtual Time Reversal Method for Passive Direction of Arrival Estimation

The paper presents a new method to estimate the two-dimensional (2D) direction of arrival (DOA) (the azimuth angle and the elevation angle) of electromagnetic signal emitted from a single communication station. This method is passive and accurate in the case of low signal-noise ratio (SNR) based on the virtual time reversal (VTR) theory. In order to illustrate its principle, the theoretical formulas of VTR direction finding with uniform circular array (UCA) are derived firstly. Based on these formulas, the implementation scheme for estimating azimuth angle and elevation angle passively is then provided. In the derivation, the strict mathematical proof for compressing planar search area to a curve line is proposed, reducing the complexity of VTR algorithm greatly. Finally, the simulation experiments are performed to validate the performance of VTR algorithm.The results show that the VTR method is effective and it delivers accurate DOA estimation in the case of low SNR.


Introduction
It is very important to detect the direction of concealed enemy communication station passively and accurately under the condition of electronic warfare (EW).Until now, several commonly used methods for estimating the direction of arrival (DOA) have been developed, such as phase comparison direction finding method [1], multiple signal classification (MUSIC) [2], and estimation of signal parameters via rotational invariance techniques (ESPRIT) [3].The performance of these methods in conventional electromagnetic environment is satisfactory, but it is not quite enough to estimate the DOA of weak electromagnetic signal in the case of modern EW environment.For example, the performance of phase comparison method is often affected by low signalnoise radio (SNR) levels, and there exists the ambiguity of bearing.The MUSIC method has high resolution, but it is difficult to implement in real time under the case of large sample.In order to improve the computation complexity, a hybrid MUSIC approach with uniform circular array (UCA) is proposed in [4].The results show that the average runtime of this combined algorithm is reduced, but the error of elevation angle is almost 5 ∘ when the SNR is 0 dB, which is not satisfactory.A sparse MUSIC algorithm with UCA is provided in [5], and it gives more accurate DOA estimation than traditional MUSIC.However, it is impossible to get a root-mean-square error (RMSE) below 1 ∘ at SNR = 0 dB when the number of elements is less than 10.ESPRIT offers advantages over MUSIC algorithm by avoiding the spectrum peak search and reducing the effect of sensor performance variability.However, it cannot be used for UCA directly, and the estimation error of ESPRIT is greater than MUSIC under normal circumstances [6,7].Therefore developing a new method to estimate the DOA of weak electromagnetic signal accurately is the necessary to EW.
The time reversal (TR) direction finding method has some advantages, such as simple algorithm, accurate estimation result, and self-focusing without prior knowledge about medium characteristic and array distribution [8], so it has a potential ability to estimate the DOA of weak electromagnetic signal.It is worth mentioning that the TR technology accomplishes self-refocusing without knowing the distribution of TR array; therefore it is no limitation in the geometry of TR array.TR theory [9] is proposed by Fink in 1989, and it has been developed in various areas, such as nondestructive testing [10], lithotripsy [11], acoustics [12], imaging [13], civil engineering [14], radar [15], communications [16], and electromagnetic [17][18][19].Although the focusing characteristic of TR method has been studied extensively, most of published papers applied the classic idea proposed by Fink [20].In the classical TR theory, a physical backpropagation process for time-reversed signal is necessary, and a transducers array which is placed at the source to receive the backpropagated signal is needed for observing the refocusing process.Hence it is unsuitable to detect the direction of the enemy communication station because we cannot set up a receive array in the position of enemy communication station, and the physical backpropagation process for time-reversed signal is meaningless.
In order to solve the problem above, a virtual time reversal (VTR) passive direction finding method used for single communication station is proposed, which works with TR array only, and without receive array at source, and the physical backpropagation process for time-reversed signal is replaced by the virtual calculation.For verifying the performance of the proposed method, some experiments are performed.The simulation results show that the VTR DOA estimation method has a smaller root-mean-square error (RMSE) at low SNR levels than phase comparison method, spectral MUSIC, and real beam space MUSIC (RB-MUSIC).
The remainder of this paper is organized as follows.Section 2 proposes the VTR passive direction finding method for estimating azimuth angle and elevation angle.In Section 3, the time delay correction is given, and the steps of VTR algorithm are summarized.The simulation experiment results and performance analyses of the VTR algorithm are shown in Section 4. We summarize our works in Section 5.

Fundamental of VTR Direction Finding Method
In this section, the direction finding method of azimuth angle is firstly given, and then under known azimuth angle the direction finding method of elevation angle is proposed.

Model and Algorithm.
With no limitation in geometry of TR array, we can establish the coordinate system as shown in Figure 1.The TR array is an antenna array, and it only receives signal.The antenna elements, assumed to be identical and omnidirectional, are uniformly distributed over the circumference of a circle of radius  0 , and the number of antennas is 8.A cylindrical coordinate system is used to represent the arrival direction of incoming electromagnetic wave.The origin of the coordinate system is located at the centre of the array.The position of every element is  0 ( 0 ,  0 ,  0 ), where  0 = 20 m,  0 = ( − 1) × 45 ∘ ,  = 1, 2, . . ., 8,  0 = 0. Suppose that the coordinate of communication station is (  ,   ,   ), and the projection in  plane is   (  ,   , 0), as shown in Figure 1.Obviously the azimuth angle of projection is the same as source.The distance from source to th antenna element is defined as    ; according to cosine theorem we have Hence, the time delay from  to th element is given by where  = 3 × 10 8 m/s is the velocity of the electromagnetic wave.
Since wave propagation in electromagnetic environment only changes time shift and amplitude reduction of initial source wave, the received signal of th element emanating from source by () is given by The duration of each signal recorded by transducer is .After time-reversed the received signal of th element, the output signal is written as where  ∈ [0, ].
We retransmit time-reversed signals and compensate the forward transmission attenuation of () via virtual computation; then the retransmitted signals propagate back and do superposition at source.The received signal at source can be expressed by So, the energy of received signal at source is given by According to the analysis above, the energy () is maximized at source.In other words, the TR algorithm accomplishes refocusing on the source.
In fact, the refocusing of signals is conditional, and we need to know the amplitude attenuation factor 1/   in advance.However, the position of enemy communication station is unknown making it difficult to obtain the amplitude attenuation factor.In order to solve this problem, we set channel fading gently; thus the amplitude attenuation factor can be considered a constant (it is not changing with time) in the record time .Hence, compensating or not does not influence the peak point of ().In view of the reasons above, we can normalize the signal energy in VTR process, without considering the amplitude attenuation compensation.

VTR Process for Estimating Azimuth Angle.
In EW, it is impossible to place a receive array at enemy communication station.So, we must use virtual computation to retransmit and backpropagate the time-reversed signal instead of emitting time-reversed signal to medium truly.The VTR process is as follows: antenna array receives and records signal emitted by enemy communication station firstly, and then this signal is time reversed.After that we retransmit the time-reversed signal virtually via different time delays and conclude the virtual signal energy of all scanning points in the given search area and then find the maximum value point of refocused signal energy.And the azimuth angle of this point is the azimuth angle of source.In order to expound the principle of VTR process,  plane is set as the search area, as shown in Figure 2. It shows that the search area is partitioned by concentric circles and rays from centre of UCA.The distance between adjacent two concentric circles is Δ, and the radius of each concentric circle is   = Δ ( = 1,2,...).The included angle of adjacent two rays is Δ, and the azimuth angle of each ray is   = Δ ( = 0, 1, . . ., (360 ∘ /Δ) − 1).
The intersection points of concentric circles and rays are the scanning points.Suppose the coordinate of th scanning point is   (  ,   , 0).According to the established process of (1), we have the distance    from this scanning point to th element: where  = 1, 2, . . ., 8.
The time delay from th scanning point to th element is expressed by Without considering the amplitude attenuation of transmission, the received signal of th element is given by Time reversed the received signal in record time , and we get the time-reversed signal as follows: where  ∈ [0, ].
We retransmit the time-reversed signal virtually via time delays for different scanning points, and every channel signal does superposition at   ; then we get the virtual received signal at   : After above, the virtual signal energy at   is given by Obviously,   () is the function of   , and different scanning points have different energy values.We can attest that the maximum energy value is only relevant to the azimuth angle of scanning point under certain conditions, and this azimuth angle equals the azimuth angle of projection of source in  plane.Based on this conclusion we can simplify the planar search area which is relevant to  and  to a circular line which is only relevant to , as shown in Figure 3.The mathematical proof process of this conclusion is seen as in Figure 3.
Suppose the form of source signal is sine; taking advantage of trigonometric function, (11) is simplified as P il (r i ,  l , 0) P il (r i ,  l , 0) r i r i S  (r s ,  s , 0) S  (r s ,  s , 0) Let  =  2 0 +  2  ,  = 2 0   ,  =  2 0 +  2  ,  = 2 0   , and according to (1), ( 7), we have ) . ( In order to simplify (14), we can expand the radical sign term by Taylor expansion, which needs |(± cos(( − 1) ⋅ 45 ∘ −   ) +   14) is expressed by Similarly, we have According to ( 15)-( 16), (13) becomes where Obviously, the maximum value of   () is consistent with   (), and we can get the maximum value of   () by means of solving derivative of   () with respect to   .Because the sine function term has nothing to do with   , solving derivative of   () becomes solving derivative of  1 .
Let  = − 1 cos   +  1 cos   ,  = − 1 sin   +  1 sin   .Equation (18) becomes We have known And then (20) can be given by Solving the derivative of  1 with respect to   and letting  1 /  = 0, we can get cos   sin   = cos   sin   . ( From the above mathematical proof, we can see that only if   =   ,   () is maximized.Hence, the correctness of VTR passive direction finding method about azimuth angle can be confirmed under the certain condition of   < /4.But when   ≥ /4, the Taylor expansion of the radical sign term is unavailable, and the proposed method is also unavailable.In addition, the maximum value of   () is only relevant with   and has nothing to do with what   is when   < /4.So the azimuth angle obtained by scanning along with a circular line whose radius is   ( 0 ≪   ) is the same as the azimuth angle obtained by scanning the whole planar area.Thus the number of scanning points is reduced greatly, and time consumption of the proposed VTR method is also reduced greatly.
The azimuth angle is obtained by scanning a circular line.According to this azimuth angle the following section develops the VTR process for estimating the elevation angle.

VTR Process for Estimating Elevation
Angle.An implementation scheme for estimating the elevation angle with UCA as shown in Figure 1 is proposed in this section.We have obtained the azimuth angle   of enemy communication station in Section 2.2, and then we can simplify the search area of elevation angle detection to a part of    plane, as shown in Figure 4.
Figure 4 shows that   =   cos   ,   =   sin   ; then the distance    from source to th element in ( 1) is given by where  = 1, 2, . . ., 8,   < 45 ∘ .Figure 4 also tells us that the search area is a quarter of    plane.Suppose the distance between two adjacent circular arcs is Δ, and the radius of th circular arc is   = Δ (Δ = 500 m,  = 1, 2, . ..).The elevation angle of ray from centre of UCA is increased with an increment Δ = 1 ∘ , and the elevation angle of th ray is   = Δ ( = 0, 1, . . ., 89).Set the intersection points of th circular arc with th ray as scanning point   (  ,   ),   =   cos   ; then the distance from   to th element is given by Further, using the VTR process given in Section 2.2 to scan virtual sources in the whole search area of elevation angle we can obtain the elevation angle of source.However, it still involves the time-consuming problem while the search area of elevation angle is a quarter of    plane.So simplifying the search area is necessary.The mathematical proof of simplifying search area is given in what follows.
Similar to (11), we can get the superimposed signals at   as Suppose the source signal is sine form; taking advantage of the trigonometric function transformation we can simplify (25) as follows: Let  =  2 0 +  2  ,  = 2 0   ,  =  2 0 +  2  ,  = 2 0   .According to (23) and (24), we have Similarly, we have According to ( 27)-( 28), (26) becomes where Obviously, the maximum value of   () is consistent with  11 , and we can get the maximum value of   () by means of solving derivative of   () with respect to   .Because the sine function term has nothing to do with   , solving derivative of   () becomes solving derivative of  11 .

Because of 𝛼/𝛼 = 𝑟
From the above mathematical proof, we can know that when   =   ,   () is maximized, the correctness of passive elevation angle detection method based on VTR theory given by this section is confirmed.In addition, the maximum value of   () is only relevant with   and has nothing to do with what   is.So the result obtained by scanning along with any circular arc whose radius is   (  ≫  0 ) is the same as the result obtained by scanning the whole area.This conclusion simplifies the sector search area to a circular arc, as shown in Figure 5.So this conclusion reduces the number of scanning points greatly, so that the time consumption of elevation angle detecting is reduced greatly.
After the above analyses, the simplifying process about search area is finished, and the simplified search area only contains a circular line (see Figure 3) and a circular arc (see Figure 5).Compared with the search area described in Figures 2 and 4, the number of scanning points is reduced greatly, so the time consumption in VTR direction finding method is reduced greatly, and the algorithm complexity is also reduced greatly.
With the description about the fundamental of VTR direction finding method the algorithm of VTR passive direction finding method is given in the following section.

Algorithm of VTR Direction
Finding Method If the frequency of signal converts from RF to IF, the phase of RF signal also converts to IF signal coincident with the frequency conversion.Hence, the propagating time delay  of RF signal is changed.We must correct  during VTR process if we want to obtain the direction of scanning point exactly.Suppose the corrected time delay is   , and we have the following expression: Thus, where  is the phase of input signal,   is the RF of input signal,   is the IF after frequency conversion, and  =   /  is the corrective factor of time delay.From (37) we can see that  is infinity if   = 0, so the time delay correction is failure for zero IF signal.Thus VTR direction finding method is suitable for signal whose converted frequency is not zero IF.

VTR Algorithm for Estimating DOA.
We have estimated the azimuth angle and elevation angle in Sections 2.2 and 2.3 and corrected the time delay in Section 3.1.Thus we summarize the steps of VTR direction finding algorithm as follows.
Step 1.The receive antenna array receives electromagnetic wave signal from enemy communication station and records it in  seconds.
Step 2. Time reversed the recorded signal.
Step 3. Determine a circular line search area of azimuth angle and calculate the time delays from scanning points to each antenna element in search area and then translate them into sample points using (38) as follows: where   denotes translated sample points in search area of azimuth angle and   is the sample rate.

Mathematical Problems in Engineering
If we convert a RF signal to IF signal, the time delays from scanning points to each antenna element must be corrected by making use of (37); thus (38) will become where     is the corrective time delay.The time-reversed signal is retransmitted to each scanning point via the calculated sample points.And retransmitted signal by each antenna element is added together and the signal energy is calculated in each scanning point.
Step 4. Find the maximum value of all signal energy and obtain the azimuth angle of enemy communication station.
Step 5.According to the estimated azimuth angle, we determine the circular arc search area of elevation angle and calculate the time delays from scanning points to each antenna element in search area and then translate them into sample points using the following equation: where   denotes translated sample points in search area of elevation angle.
Similarly, if a RF signal is converted to a IF signal, the time delays from scanning points to each antenna element must be corrected by making use of (37); thus (40) will become where     is the corrective time delay.The time-reversed signal in Step 2 is retransmitted to each scanning point via the calculated sample points by (41).And retransmitted signal by each antenna element is added together and the signal energy is calculated in each scanning point.
Step 6. Find the maximum value of all signal energy in search area of elevation angle and obtain the elevation angle of enemy communication station.
From the above analyses of VTR direction finding method we can observe that the principle of the proposed VTR method is the same as classical TR method.The only difference is the working model for retransmitting the time-reversed signal to medium.In the classical TR method, the time-reversed signal is retransmitted to medium truly, and then, the retransmitted signal is received by the array which located at source and refocus at source.But the VTR method retransmits time-reversed signal to the given search area virtually via time delays and refocus at source and needs not the receive array located at source.So the performance of two methods is the same, but VTR method can be used to estimate the emission source passively and the classical TR method cannot.What is more important is that we simplify the search area in our VTR method, and this idea is never mentioned in classical TR method.
The following section documents the results of computer simulations of the VTR method.

Computer Simulations and Performance Analyses
Consider a UCA as shown in Figure 1, and all numerical results given in this section pertain to this array configuration.Consider SystemView simulator as the simulation tool.
In order to demonstrate and evaluate the new algorithm proposed in this paper, the specific parameters are set as follows.
(1) Parameters of source: the location of source is  (20 km, 22 ∘ , 7.28 km), and the elevation angle is   = 20 ∘ < 45 ∘ .Consider the source which is modulated by BPSK emitting 300 kbps pseudorandom bit sequences.The carrier frequency is 10.7 MHz, and the amplitude is 1 V.
(2) Parameters of channel: the signal is received in the presence of additive white Gaussian noise.The SNR level is 0 dB (SNR = 10 lg   /(( 0 /2) ⋅ BW),   is the signal power,  0 /2 is the noise power spectrum density, and BW is the receive bandwidth of signal).
And the transmission amplitude attenuation cannot be simulated.
(3) Parameters of received signal: consider the receive bandwidth which is 600 kHz; the IF of received signal is   = 700 kHz.The length of window for recording signal is  = 4 ms and the centre of window is slidable.
(5) Parameters of DOA searching: the radius of circular line is   = 20 km, and the scanning increment of azimuth angle is Δ = 1 ∘ .
The virtual received signal of scanning point is produced by (10) and (11).The energy value is calculated by (12), and Figure 6 shows the energy of scanning points.In Figure 6 the top is the energy curve of azimuth angle from 0 ∘ to 359 ∘ , the middle is the energy distribution along the circular line, and the bottom is the planform of the middle.From this figure we can see that the corresponding angle of peak value is 22 ∘ , and it is coincident with the azimuth angle of source.
In order to observe the estimation error of azimuth angle in the whole search area, we set the radius of circular line to be   = Δ (Δ = 500 m,  = 1, 2, . . ., 40), and Figure 7 shows the RMSEs for 40 scanning circles.In simulation, the energy calculation is repeated for 100 times, which is a Monte-Carlo simulation of 100 runs for different radiuses.We can see from Figure 7 that the RMSE of any radius is close to others, which reveals that the estimation error is relatively stable when the radius of circular line is changed.Figure 7 also shows that the VTR passive direction finding method has a relatively small bearing estimation error when the SNR level of input signal is 0 dB.
In the second experiment we study the estimation error of azimuth angle for different SNR levels and different record time .We consider the radius of circular line is   = 20 km, the record time is  = 4 ms, and set different power spectral  densities of additive white Gaussian noise to change the SNR levels.In Figure 8 the RMSEs obtained by proposed algorithm are compared to the RMSEs obtained by means of the phase comparison method, spectral MUSIC, and RB-MUSIC.Reference [7] also describes a UCA-ESPRIT estimator, but it has lower estimation accuracy than RB-MUSIC estimator, which is not considered here.The phase mode of RB-MUSIC is set as 3.The result of phase comparison method is the average value of two pairs of baselines (baseline 1, 5 and baseline 2, 6, baseline 3, 7 and baseline 4, 8. Using each pair of baselines, we can obtain the azimuth angle and the elevation angle).We can see from Figure 8 that the RMSEs obtained by proposed method are less than the compared methods at low SNR levels.But the estimation accuracy of proposed method is lower than spectral MUSIC and RB-MUSIC when SNR > 14 dB.This result indicates that VTR passive direction finding method delivers accurate azimuth angle estimation for weak signal.In order to study the estimation error for different record time , we consider the radius of scanning circle as   = 20km and set the SNR level of the input signal as 0 dB; Figure 9 plots the RMSEs obtained by VTR algorithm and compared methods for different record time.
From Figure 9 we can observe that the performance of DOA estimation becomes better in collaboration with  increasing, and the estimation accuracy of VTR algorithm is better than compared methods at SNR = 0 dB.
virtual time reversal approach that this approach gives us a significant result for simplifying the search area.Except that, the virtual time reversal approach uses the time delay compensation instead of the physical backpropagation of time-reversed signals, which changes the working model of the classical time reversal theory.For all shapes of received arrays the time delay compensation process can be implemented, so the virtual time reversal approach is not limited to the geometry of received array.For multiple sources, the channelized receiver, which is used to divide the received signal in frequency, is needed before DOA estimation, and then the virtual time reversal approach proposed can be used to estimate the DOAs for multiple sources.

Conclusions
This paper presents a new method to estimate the DOA of electromagnetic wave signal emanating from a single communication station passively and accurately based on the VTR theory.The advantages of this method are the simpler algorithm, the higher precision, and the concealed working model.This method can accomplish the DOA estimation of single communication station passively and accurately under the condition of   < 45 ∘ .Based on the mathematical models of VTR direction finding algorithm described in this paper, a significant conclusion for simplifying the planar search area to a curve line is given by strict mathematical proof.And this conclusion reduces the number of scanning points and time consumption in VTR direction finding method greatly.After that the formula of time delay correction when a RF signal is converted to an IF signal is derived and the steps of VTR algorithm are provided.In order to evaluate the performance of this method for estimating DOA, the simulation experiments are performed.The results show that the proposed method has better estimation accuracy for the weak signal than phase comparison method, the spectral MUSIC, and the RB-MUSIC.
It is worth noting that the antenna array used in this paper is placed in the same horizontal plane and worked passively, so it is concealed.Thus the proposed method is more suitable to detect the direction of enemy communication station passively.It also can be used to locate the communication station passively if we establish two antenna arrays and let them work together.The next challenges are to estimate the DOAs of multiple sources and coherent sources based on the VTR theory.

Figure 2 :
Figure 2: Search area of azimuth angle.

Figure 3 :
Figure 3: Simplified search area of azimuth angle.

Figure 4 :
Figure 4: Search area of elevation angle.

Figure 5 :
Figure 5: Simplified search area of elevation angle.
Solving the derivative of  11 with respect to   and setting  11 /  = 0, we can get