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

Authors: Yongqing Fu; Wei Liu; Ruijie Bai; Jingrui Li; Jinlin Wang
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

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.


---

## Body

## 1. 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 signal-noise 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–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. Section2 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.


## 2. 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.


### 2.1. Model and Algorithm

With no limitation in geometry of TR array, we can establish the coordinate system as shown in Figure1. 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 r
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 S
0
r
0
,
θ
0
k
,
z
0, where r
0
=
20 m, θ
0
k
=
k
-
1
×
4
5
∘, k
=
1,2
,
…
,
8,  z
0
=
0.Figure 1
VTR receive array.Suppose that the coordinate of communication station isS
r
s
,
θ
s
,
z
s, and the projection in r
o
θ plane is S
′
r
s
,
θ
s
,
0, as shown in Figure 1. Obviously the azimuth angle of projection is the same as source. The distance from source to kth antenna element is defined as R
S
k; according to cosine theorem we have
(1)
R
S
k
=
r
0
2
+
r
s
2
-
2
r
s
r
0
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
z
s
2
1
/
2
,
k
=
1,2
,
…
,
8
.Hence, the time delay fromS to kth element is given by
(2)
τ
S
k
=
R
S
k
c
,
where c
=
3
×
1
0
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 ofkth element emanating from source by S
(
t
) is given by
(3)
S
k
t
=
S
t
-
τ
S
k
R
S
k
.The duration of each signal recorded by transducer isT. After time-reversed the received signal of kth element, the output signal is written as
(4)
S
¯
k
t
=
S
k
T
-
t
=
S
T
-
t
-
τ
S
k
R
S
k
,
where t
∈
0
,
T.We retransmit time-reversed signals and compensate the forward transmission attenuation ofS
(
t
) via virtual computation; then the retransmitted signals propagate back and do superposition at source. The received signal at source can be expressed by
(5)
S
¯
t
=
∑
k
=
1
8
R
S
k
S
¯
k
t
-
R
S
k
c
=
∑
k
=
1
8
R
S
k
S
T
-
t
R
S
k
=
8
S
T
-
t
.So, the energy of received signal at source is given by(6)
E
t
=
S
¯
2
t
=
8
S
T
-
t
2
=
64
S
2
T
-
t
.According to the analysis above, the energyE
(
t
) 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 factor1
/
R
S
k 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 T. Hence, compensating or not does not influence the peak point of E
(
t
). In view of the reasons above, we can normalize the signal energy in VTR process, without considering the amplitude attenuation compensation.


### 2.2. 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,r
o
θ 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 Δ
r, and the radius of each concentric circle is r
i
=
i
Δ
r
(
i
=
1,2
,
…
). The included angle of adjacent two rays is Δ
θ, and the azimuth angle of each ray is θ
l
=
l
Δ
θ
(
l
=
0,1
,
…
,
(
36
0
∘
/
Δ
θ
)
-
1
). The intersection points of concentric circles and rays are the scanning points.Figure 2
Search area of azimuth angle.Suppose the coordinate ofi
lth scanning point is P
i
l
r
i
,
θ
l
,
0. According to the established process of (1), we have the distance R
P
k
i
l from this scanning point to kth element:
(7)
R
P
k
i
l
=
r
0
2
+
r
i
2
-
2
r
i
r
0
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
1
/
2
,
where k
=
1,2
,
…
,
8.The time delay fromi
lth scanning point to kth element is expressed by
(8)
τ
P
k
i
l
=
R
P
k
i
l
c
.Without considering the amplitude attenuation of transmission, the received signal ofkth element is given by
(9)
S
k
t
=
S
t
-
τ
S
k
.Time reversed the received signal in record timeT, and we get the time-reversed signal as follows:
(10)
S
¯
k
t
=
S
k
T
-
t
=
S
T
-
t
-
τ
S
k
,
where t
∈
0
,
T.We retransmit the time-reversed signal virtually via time delays for different scanning points, and every channel signal does superposition atP
i
l; then we get the virtual received signal at P
i
l:
(11)
S
¯
i
l
t
=
∑
k
=
1
8
S
¯
k
t
-
τ
P
k
i
l
=
∑
k
=
1
8
S
T
-
t
+
τ
P
k
i
l
-
τ
S
k
.After above, the virtual signal energy atP
i
l is given by
(12)
E
i
l
t
=
S
¯
i
l
t
2
=
∑
k
=
1
8
S
¯
k
t
-
τ
P
k
i
l
2
=
∑
k
=
1
8
S
T
-
t
+
τ
P
k
i
l
-
τ
S
k
2
.Obviously,E
i
l
(
t
) is the function of P
i
l, 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 r
o
θ plane. Based on this conclusion we can simplify the planar search area which is relevant to r 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.Figure 3
Simplified search area of azimuth angle.Suppose the form of source signal is sine; taking advantage of trigonometric function, (11) is simplified as
(13)
S
¯
i
l
t
=
2
∑
k
=
1
4
sin
T
-
t
+
R
P
k
i
l
+
R
P
k
+
4
i
l
-
R
S
k
+
R
S
k
+
4
2
c
aaaaaaaaa
·
cos
⁡
R
P
k
i
l
-
R
P
k
+
4
i
l
-
R
S
k
-
R
S
k
+
4
2
c
.Leta
=
r
0
2
+
r
i
2, b
=
2
r
0
r
i, a
-
=
r
0
2
+
r
s
2, b
-
=
2
r
0
r
s, and according to (1), (7), we have
(14)
R
S
k
+
R
S
k
+
4
=
a
-
-
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
z
s
2
1
/
2
+
a
-
+
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
z
s
2
1
/
2
=
a
-
1
+
-
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
z
s
2
a
-
1
/
2
aaaaa
+
1
+
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
z
s
2
a
-
1
/
2
.In order to simplify (14), we can expand the radical sign term by Taylor expansion, which needs ±
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
z
s
2
/
a
-
<
1. So, we can get b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
a
-
+
z
s
2
/
a
-
<
1, and the demand is met. Because b
-
/
a
-
=
2
r
0
/
r
s
+
r
0
2
/
r
s
≪
1, and 0
≤
z
s
2
/
a
-
<
1, we have 0
≤
r
s
tan
φ
s
2
/
r
0
2
+
r
s
2
<
1, where φ
s is the elevation angle of source, and φ
s
∈
0
,
π
/
2.From0
≤
ta
n
2
φ
s
/
1
+
r
0
2
/
r
s
2
<
1, we have 0
≤
tan
φ
s
<
1; that is, 0
≤
φ
s
<
π
/
4. So, when 0
≤
φ
s
<
π
/
4, (14) is expressed by
(15)
R
S
k
+
R
S
k
+
4
≈
a
-
1
-
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
a
-
-
z
s
2
/
a
-
2
a
a
a
a
a
a
a
a
+
1
+
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
a
-
+
z
s
2
/
a
-
2
=
2
a
-
+
z
s
2
a
-
.Similarly, we have(16)
R
P
k
i
l
+
R
P
k
+
4
i
l
=
a
-
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
1
/
2
+
a
+
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
1
/
2
=
a
1
-
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
a
1
/
2
+
1
+
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
a
1
/
2
≈
a
1
-
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
2
a
+
1
+
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
2
a
=
2
a
,
R
S
k
-
R
S
k
+
4
=
-
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
a
-
,
R
P
k
i
l
-
R
P
k
+
4
i
l
=
-
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
a
.According to (15)-(16), (13) becomes
(17)
S
¯
i
l
t
=
2
sin
T
-
t
+
2
a
-
2
a
-
+
z
s
2
/
a
-
2
c
S
¯
1
,
where
(18)
S
¯
1
=
∑
k
=
1
4
cos
⁡
-
ρ
1
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
a
a
a
a
a
a
a
a
a
a
a
+
ρ
-
1
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
(19)
ρ
1
=
b
2
c
a
=
r
0
c
1
+
r
0
/
r
i
2
1
/
2
≈
r
0
1
-
r
0
/
r
i
2
/
2
c
ρ
-
1
=
b
-
2
c
a
-
=
r
0
c
1
+
r
0
/
r
s
2
1
/
2
≈
r
0
1
-
r
0
/
r
s
2
/
2
c
.Obviously, the maximum value ofE
i
l
(
t
) is consistent with S
¯
i
l
t, and we can get the maximum value of S
¯
i
l
t by means of solving derivative of S
¯
i
l
t with respect to θ
l. Because the sine function term has nothing to do with θ
l, solving derivative of S
¯
i
l
t becomes solving derivative of S
¯
1.LetX
=
-
ρ
1
cos
⁡
θ
l
+
ρ
-
1
cos
⁡
θ
s, Y
=
-
ρ
1
sin
θ
l
+
ρ
-
1
sin
θ
s. Equation (18) becomes
(20)
S
¯
1
=
cos
⁡
X
+
cos
⁡
2
X
+
Y
2
+
cos
⁡
Y
+
cos
⁡
2
X
-
Y
2
.We have knownr
0
≪
r
i, r
0
≪
r
s, ρ
1
≪
1, ρ
-
1
≪
1, so X
≪
1, Y
≪
1. And then (20) can be given by
(21)
S
¯
1
=
1
-
X
2
2
+
1
-
2
X
+
Y
2
/
2
2
+
1
-
Y
2
2
+
1
-
2
X
-
Y
2
/
2
2
=
4
-
X
2
-
Y
2
.Solving the derivative ofS
¯
1 with respect to θ
l and letting d
S
¯
1
/
d
θ
l
=
0, we can get
(22)
cos
⁡
θ
s
sin
θ
l
=
cos
⁡
θ
l
sin
θ
s
.From the above mathematical proof, we can see that only ifθ
l
=
θ
s, S
¯
i
l
t is maximized. Hence, the correctness of VTR passive direction finding method about azimuth angle can be confirmed under the certain condition of φ
s
<
π
/
4. But when φ
s
≥
π
/
4, the Taylor expansion of the radical sign term is unavailable, and the proposed method is also unavailable. In addition, the maximum value of S
¯
i
l
t is only relevant with θ
l and has nothing to do with what r
i is when φ
s
<
π
/
4. So the azimuth angle obtained by scanning along with a circular line whose radius is r
i
(
r
0
≪
r
i
) 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.


### 2.3. VTR Process for Estimating Elevation Angle

An implementation scheme for estimating the elevation angle with UCA as shown in Figure1 is proposed in this section. We have obtained the azimuth angle θ
s of enemy communication station in Section 2.2, and then we can simplify the search area of elevation angle detection to a part of θ
s
O
Z plane, as shown in Figure 4.Figure 4
Search area of elevation angle.Figure4 shows that r
s
=
R
s
cos
⁡
φ
s, z
s
=
R
s
sin
φ
s; then the distance R
S
k from source to kth element in (1) is given by
(23)
R
S
k
=
r
0
2
+
R
s
cos
⁡
φ
s
2
-
2
r
0
R
s
cos
⁡
φ
s
a
a
a
·
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
R
s
sin
φ
s
2
1
/
2
=
r
0
2
+
R
s
2
-
2
r
0
R
s
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
1
/
2
,
where k
=
1,2
,
…
,
8, φ
s
<
4
5
∘.Figure4 also tells us that the search area is a quarter of θ
s
O
Z plane. Suppose the distance between two adjacent circular arcs is Δ
l, and the radius of jth circular arc is R
j
=
j
Δ
l (Δ
l
=
500 m, j
=
1,2
,
…). The elevation angle of ray from centre of UCA is increased with an increment Δ
φ
=
1
∘, and the elevation angle of mth ray is φ
m
=
m
Δ
φ
(
m
=
0,1
,
…
,
89
). Set the intersection points of jth circular arc with mth ray as scanning point P
j
m
r
j
,
φ
m, r
j
=
R
j
cos
⁡
φ
m; then the distance from P
j
m to kth element is given by
(24)
R
P
k
j
m
=
r
0
2
+
R
j
cos
⁡
φ
m
2
-
2
r
0
R
j
cos
⁡
φ
m
a
a
a
·
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
R
j
sin
φ
m
2
1
/
2
=
r
0
2
+
R
j
2
-
2
r
0
R
j
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
1
/
2
.Further, using the VTR process given in Section2.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 θ
s
O
Z 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 P
j
m as
(25)
S
¯
j
m
t
=
∑
k
=
1
8
S
T
-
t
+
R
P
k
j
m
-
R
S
k
c
.Suppose the source signal is sine form; taking advantage of the trigonometric function transformation we can simplify (25) as follows:
(26)
S
¯
j
m
=
2
∑
k
=
1
4
sin
T
-
t
+
R
P
k
j
m
+
R
P
k
+
4
j
m
-
R
S
k
+
R
S
k
+
4
2
c
·
cos
⁡
R
P
k
j
m
-
R
P
k
+
4
j
m
-
R
S
k
-
R
S
k
+
4
2
c
.LetA
=
r
0
2
+
R
j
2, B
=
2
r
0
R
j, A
-
=
r
0
2
+
R
s
2, B
-
=
2
r
0
R
s. According to (23) and (24), we have
(27)
R
S
k
+
R
S
k
+
4
=
A
-
-
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
1
/
2
+
A
-
+
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
1
/
2
=
A
-
1
-
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
A
-
1
/
2
a
a
a
a
a
a
a
a
a
+
1
+
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
A
-
1
/
2
≈
A
-
1
-
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
A
-
2
a
a
a
a
a
a
.
a
a
+
1
+
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
A
-
2
=
2
A
-
.Similarly, we have(28)
R
P
k
j
m
+
R
P
k
+
4
j
m
=
A
-
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
1
/
2
+
A
+
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
1
/
2
=
A
1
-
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
A
1
/
2
+
1
+
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
A
1
/
2
≈
A
1
-
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
A
2
+
1
+
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
A
2
=
2
A
,
R
S
k
-
R
S
k
+
4
=
-
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
A
-
,
R
P
k
j
m
-
R
P
k
+
4
j
m
=
-
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
A
.According to (27)-(28), (26) becomes
(29)
S
¯
j
m
t
=
2
sin
T
-
t
+
A
-
A
-
c
S
¯
11
,
where
(30)
S
¯
11
=
∑
k
=
1
4
cos
⁡
-
α
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
a
a
a
a
a
a
+
α
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
,
(31)
α
=
B
2
c
A
=
r
0
c
1
+
r
0
/
R
j
2
1
/
2
≈
r
0
1
-
r
0
/
R
j
2
/
2
c
,
α
-
=
B
-
2
c
A
-
=
r
0
c
1
+
r
0
/
R
s
2
1
/
2
≈
r
0
1
-
r
0
/
R
s
2
/
2
c
.Obviously, the maximum value ofS
¯
j
m
t is consistent with S
¯
11, and we can get the maximum value of S
¯
j
m
t by means of solving derivative of S
¯
j
m
t with respect to φ
m. Because the sine function term has nothing to do with φ
m, solving derivative of S
¯
j
m
t becomes solving derivative of S
¯
11.LetX
=
-
α
cos
⁡
θ
s
cos
⁡
φ
m
+
α
-
cos
⁡
θ
s
cos
⁡
φ
s, Y
=
-
α
sin
θ
s
cos
⁡
φ
m
+
α
-
sin
θ
s
cos
⁡
φ
s. Equation (30) becomes
(32)
S
¯
11
=
cos
⁡
X
+
cos
⁡
2
X
+
Y
2
+
cos
⁡
Y
+
cos
⁡
2
X
-
Y
2
.We have knownr
0
≪
R
j, r
0
≪
R
s, α
≪
1, α
-
≪
1, so X
≪
1, Y
≪
1. And then (32) can be given by
(33)
S
¯
11
=
1
-
X
2
2
+
1
-
2
X
+
Y
2
/
2
2
+
1
-
Y
2
2
+
1
-
2
X
-
Y
2
/
2
2
=
4
-
X
2
-
Y
2
.Solving the derivative ofS
¯
11 with respect to φ
m and setting d
S
¯
11
/
d
φ
m
=
0, we can get
(34)
cos
⁡
φ
m
=
α
-
cos
⁡
φ
s
α
.Because ofα
-
/
α
=
r
0
1
-
r
0
/
R
s
2
/
2
/
c
/
r
0
1
-
r
0
/
R
j
2
/
2
/
c
≈
1, we have
(35)
cos
⁡
φ
m
=
cos
⁡
φ
s
.From the above mathematical proof, we can know that whenφ
m
=
φ
s, S
¯
j
m
t 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 S
¯
j
m
t is only relevant with φ
m and has nothing to do with what R
j is. So the result obtained by scanning along with any circular arc whose radius is R
j
(
R
j
≫
r
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.Figure 5
Simplified search area of elevation angle.After the above analyses, the simplifying process about search area is finished, and the simplified search area only contains a circular line (see Figure3) 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.


## 2.1. Model and Algorithm

With no limitation in geometry of TR array, we can establish the coordinate system as shown in Figure1. 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 r
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 S
0
r
0
,
θ
0
k
,
z
0, where r
0
=
20 m, θ
0
k
=
k
-
1
×
4
5
∘, k
=
1,2
,
…
,
8,  z
0
=
0.Figure 1
VTR receive array.Suppose that the coordinate of communication station isS
r
s
,
θ
s
,
z
s, and the projection in r
o
θ plane is S
′
r
s
,
θ
s
,
0, as shown in Figure 1. Obviously the azimuth angle of projection is the same as source. The distance from source to kth antenna element is defined as R
S
k; according to cosine theorem we have
(1)
R
S
k
=
r
0
2
+
r
s
2
-
2
r
s
r
0
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
z
s
2
1
/
2
,
k
=
1,2
,
…
,
8
.Hence, the time delay fromS to kth element is given by
(2)
τ
S
k
=
R
S
k
c
,
where c
=
3
×
1
0
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 ofkth element emanating from source by S
(
t
) is given by
(3)
S
k
t
=
S
t
-
τ
S
k
R
S
k
.The duration of each signal recorded by transducer isT. After time-reversed the received signal of kth element, the output signal is written as
(4)
S
¯
k
t
=
S
k
T
-
t
=
S
T
-
t
-
τ
S
k
R
S
k
,
where t
∈
0
,
T.We retransmit time-reversed signals and compensate the forward transmission attenuation ofS
(
t
) via virtual computation; then the retransmitted signals propagate back and do superposition at source. The received signal at source can be expressed by
(5)
S
¯
t
=
∑
k
=
1
8
R
S
k
S
¯
k
t
-
R
S
k
c
=
∑
k
=
1
8
R
S
k
S
T
-
t
R
S
k
=
8
S
T
-
t
.So, the energy of received signal at source is given by(6)
E
t
=
S
¯
2
t
=
8
S
T
-
t
2
=
64
S
2
T
-
t
.According to the analysis above, the energyE
(
t
) 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 factor1
/
R
S
k 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 T. Hence, compensating or not does not influence the peak point of E
(
t
). In view of the reasons above, we can normalize the signal energy in VTR process, without considering the amplitude attenuation compensation.


## 2.2. 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,r
o
θ 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 Δ
r, and the radius of each concentric circle is r
i
=
i
Δ
r
(
i
=
1,2
,
…
). The included angle of adjacent two rays is Δ
θ, and the azimuth angle of each ray is θ
l
=
l
Δ
θ
(
l
=
0,1
,
…
,
(
36
0
∘
/
Δ
θ
)
-
1
). The intersection points of concentric circles and rays are the scanning points.Figure 2
Search area of azimuth angle.Suppose the coordinate ofi
lth scanning point is P
i
l
r
i
,
θ
l
,
0. According to the established process of (1), we have the distance R
P
k
i
l from this scanning point to kth element:
(7)
R
P
k
i
l
=
r
0
2
+
r
i
2
-
2
r
i
r
0
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
1
/
2
,
where k
=
1,2
,
…
,
8.The time delay fromi
lth scanning point to kth element is expressed by
(8)
τ
P
k
i
l
=
R
P
k
i
l
c
.Without considering the amplitude attenuation of transmission, the received signal ofkth element is given by
(9)
S
k
t
=
S
t
-
τ
S
k
.Time reversed the received signal in record timeT, and we get the time-reversed signal as follows:
(10)
S
¯
k
t
=
S
k
T
-
t
=
S
T
-
t
-
τ
S
k
,
where t
∈
0
,
T.We retransmit the time-reversed signal virtually via time delays for different scanning points, and every channel signal does superposition atP
i
l; then we get the virtual received signal at P
i
l:
(11)
S
¯
i
l
t
=
∑
k
=
1
8
S
¯
k
t
-
τ
P
k
i
l
=
∑
k
=
1
8
S
T
-
t
+
τ
P
k
i
l
-
τ
S
k
.After above, the virtual signal energy atP
i
l is given by
(12)
E
i
l
t
=
S
¯
i
l
t
2
=
∑
k
=
1
8
S
¯
k
t
-
τ
P
k
i
l
2
=
∑
k
=
1
8
S
T
-
t
+
τ
P
k
i
l
-
τ
S
k
2
.Obviously,E
i
l
(
t
) is the function of P
i
l, 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 r
o
θ plane. Based on this conclusion we can simplify the planar search area which is relevant to r 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.Figure 3
Simplified search area of azimuth angle.Suppose the form of source signal is sine; taking advantage of trigonometric function, (11) is simplified as
(13)
S
¯
i
l
t
=
2
∑
k
=
1
4
sin
T
-
t
+
R
P
k
i
l
+
R
P
k
+
4
i
l
-
R
S
k
+
R
S
k
+
4
2
c
aaaaaaaaa
·
cos
⁡
R
P
k
i
l
-
R
P
k
+
4
i
l
-
R
S
k
-
R
S
k
+
4
2
c
.Leta
=
r
0
2
+
r
i
2, b
=
2
r
0
r
i, a
-
=
r
0
2
+
r
s
2, b
-
=
2
r
0
r
s, and according to (1), (7), we have
(14)
R
S
k
+
R
S
k
+
4
=
a
-
-
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
z
s
2
1
/
2
+
a
-
+
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
z
s
2
1
/
2
=
a
-
1
+
-
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
z
s
2
a
-
1
/
2
aaaaa
+
1
+
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
z
s
2
a
-
1
/
2
.In order to simplify (14), we can expand the radical sign term by Taylor expansion, which needs ±
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
z
s
2
/
a
-
<
1. So, we can get b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
a
-
+
z
s
2
/
a
-
<
1, and the demand is met. Because b
-
/
a
-
=
2
r
0
/
r
s
+
r
0
2
/
r
s
≪
1, and 0
≤
z
s
2
/
a
-
<
1, we have 0
≤
r
s
tan
φ
s
2
/
r
0
2
+
r
s
2
<
1, where φ
s is the elevation angle of source, and φ
s
∈
0
,
π
/
2.From0
≤
ta
n
2
φ
s
/
1
+
r
0
2
/
r
s
2
<
1, we have 0
≤
tan
φ
s
<
1; that is, 0
≤
φ
s
<
π
/
4. So, when 0
≤
φ
s
<
π
/
4, (14) is expressed by
(15)
R
S
k
+
R
S
k
+
4
≈
a
-
1
-
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
a
-
-
z
s
2
/
a
-
2
a
a
a
a
a
a
a
a
+
1
+
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
a
-
+
z
s
2
/
a
-
2
=
2
a
-
+
z
s
2
a
-
.Similarly, we have(16)
R
P
k
i
l
+
R
P
k
+
4
i
l
=
a
-
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
1
/
2
+
a
+
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
1
/
2
=
a
1
-
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
a
1
/
2
+
1
+
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
a
1
/
2
≈
a
1
-
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
2
a
+
1
+
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
2
a
=
2
a
,
R
S
k
-
R
S
k
+
4
=
-
b
-
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
a
-
,
R
P
k
i
l
-
R
P
k
+
4
i
l
=
-
b
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
a
.According to (15)-(16), (13) becomes
(17)
S
¯
i
l
t
=
2
sin
T
-
t
+
2
a
-
2
a
-
+
z
s
2
/
a
-
2
c
S
¯
1
,
where
(18)
S
¯
1
=
∑
k
=
1
4
cos
⁡
-
ρ
1
cos
⁡
k
-
1
·
4
5
∘
-
θ
l
a
a
a
a
a
a
a
a
a
a
a
+
ρ
-
1
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
(19)
ρ
1
=
b
2
c
a
=
r
0
c
1
+
r
0
/
r
i
2
1
/
2
≈
r
0
1
-
r
0
/
r
i
2
/
2
c
ρ
-
1
=
b
-
2
c
a
-
=
r
0
c
1
+
r
0
/
r
s
2
1
/
2
≈
r
0
1
-
r
0
/
r
s
2
/
2
c
.Obviously, the maximum value ofE
i
l
(
t
) is consistent with S
¯
i
l
t, and we can get the maximum value of S
¯
i
l
t by means of solving derivative of S
¯
i
l
t with respect to θ
l. Because the sine function term has nothing to do with θ
l, solving derivative of S
¯
i
l
t becomes solving derivative of S
¯
1.LetX
=
-
ρ
1
cos
⁡
θ
l
+
ρ
-
1
cos
⁡
θ
s, Y
=
-
ρ
1
sin
θ
l
+
ρ
-
1
sin
θ
s. Equation (18) becomes
(20)
S
¯
1
=
cos
⁡
X
+
cos
⁡
2
X
+
Y
2
+
cos
⁡
Y
+
cos
⁡
2
X
-
Y
2
.We have knownr
0
≪
r
i, r
0
≪
r
s, ρ
1
≪
1, ρ
-
1
≪
1, so X
≪
1, Y
≪
1. And then (20) can be given by
(21)
S
¯
1
=
1
-
X
2
2
+
1
-
2
X
+
Y
2
/
2
2
+
1
-
Y
2
2
+
1
-
2
X
-
Y
2
/
2
2
=
4
-
X
2
-
Y
2
.Solving the derivative ofS
¯
1 with respect to θ
l and letting d
S
¯
1
/
d
θ
l
=
0, we can get
(22)
cos
⁡
θ
s
sin
θ
l
=
cos
⁡
θ
l
sin
θ
s
.From the above mathematical proof, we can see that only ifθ
l
=
θ
s, S
¯
i
l
t is maximized. Hence, the correctness of VTR passive direction finding method about azimuth angle can be confirmed under the certain condition of φ
s
<
π
/
4. But when φ
s
≥
π
/
4, the Taylor expansion of the radical sign term is unavailable, and the proposed method is also unavailable. In addition, the maximum value of S
¯
i
l
t is only relevant with θ
l and has nothing to do with what r
i is when φ
s
<
π
/
4. So the azimuth angle obtained by scanning along with a circular line whose radius is r
i
(
r
0
≪
r
i
) 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.


## 2.3. VTR Process for Estimating Elevation Angle

An implementation scheme for estimating the elevation angle with UCA as shown in Figure1 is proposed in this section. We have obtained the azimuth angle θ
s of enemy communication station in Section 2.2, and then we can simplify the search area of elevation angle detection to a part of θ
s
O
Z plane, as shown in Figure 4.Figure 4
Search area of elevation angle.Figure4 shows that r
s
=
R
s
cos
⁡
φ
s, z
s
=
R
s
sin
φ
s; then the distance R
S
k from source to kth element in (1) is given by
(23)
R
S
k
=
r
0
2
+
R
s
cos
⁡
φ
s
2
-
2
r
0
R
s
cos
⁡
φ
s
a
a
a
·
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
R
s
sin
φ
s
2
1
/
2
=
r
0
2
+
R
s
2
-
2
r
0
R
s
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
1
/
2
,
where k
=
1,2
,
…
,
8, φ
s
<
4
5
∘.Figure4 also tells us that the search area is a quarter of θ
s
O
Z plane. Suppose the distance between two adjacent circular arcs is Δ
l, and the radius of jth circular arc is R
j
=
j
Δ
l (Δ
l
=
500 m, j
=
1,2
,
…). The elevation angle of ray from centre of UCA is increased with an increment Δ
φ
=
1
∘, and the elevation angle of mth ray is φ
m
=
m
Δ
φ
(
m
=
0,1
,
…
,
89
). Set the intersection points of jth circular arc with mth ray as scanning point P
j
m
r
j
,
φ
m, r
j
=
R
j
cos
⁡
φ
m; then the distance from P
j
m to kth element is given by
(24)
R
P
k
j
m
=
r
0
2
+
R
j
cos
⁡
φ
m
2
-
2
r
0
R
j
cos
⁡
φ
m
a
a
a
·
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
+
R
j
sin
φ
m
2
1
/
2
=
r
0
2
+
R
j
2
-
2
r
0
R
j
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
1
/
2
.Further, using the VTR process given in Section2.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 θ
s
O
Z 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 P
j
m as
(25)
S
¯
j
m
t
=
∑
k
=
1
8
S
T
-
t
+
R
P
k
j
m
-
R
S
k
c
.Suppose the source signal is sine form; taking advantage of the trigonometric function transformation we can simplify (25) as follows:
(26)
S
¯
j
m
=
2
∑
k
=
1
4
sin
T
-
t
+
R
P
k
j
m
+
R
P
k
+
4
j
m
-
R
S
k
+
R
S
k
+
4
2
c
·
cos
⁡
R
P
k
j
m
-
R
P
k
+
4
j
m
-
R
S
k
-
R
S
k
+
4
2
c
.LetA
=
r
0
2
+
R
j
2, B
=
2
r
0
R
j, A
-
=
r
0
2
+
R
s
2, B
-
=
2
r
0
R
s. According to (23) and (24), we have
(27)
R
S
k
+
R
S
k
+
4
=
A
-
-
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
1
/
2
+
A
-
+
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
1
/
2
=
A
-
1
-
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
A
-
1
/
2
a
a
a
a
a
a
a
a
a
+
1
+
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
A
-
1
/
2
≈
A
-
1
-
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
A
-
2
a
a
a
a
a
a
.
a
a
+
1
+
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
A
-
2
=
2
A
-
.Similarly, we have(28)
R
P
k
j
m
+
R
P
k
+
4
j
m
=
A
-
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
1
/
2
+
A
+
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
1
/
2
=
A
1
-
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
A
1
/
2
+
1
+
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
A
1
/
2
≈
A
1
-
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
A
2
+
1
+
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
/
A
2
=
2
A
,
R
S
k
-
R
S
k
+
4
=
-
B
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
A
-
,
R
P
k
j
m
-
R
P
k
+
4
j
m
=
-
B
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
A
.According to (27)-(28), (26) becomes
(29)
S
¯
j
m
t
=
2
sin
T
-
t
+
A
-
A
-
c
S
¯
11
,
where
(30)
S
¯
11
=
∑
k
=
1
4
cos
⁡
-
α
cos
⁡
φ
m
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
a
a
a
a
a
a
+
α
-
cos
⁡
φ
s
cos
⁡
k
-
1
·
4
5
∘
-
θ
s
,
(31)
α
=
B
2
c
A
=
r
0
c
1
+
r
0
/
R
j
2
1
/
2
≈
r
0
1
-
r
0
/
R
j
2
/
2
c
,
α
-
=
B
-
2
c
A
-
=
r
0
c
1
+
r
0
/
R
s
2
1
/
2
≈
r
0
1
-
r
0
/
R
s
2
/
2
c
.Obviously, the maximum value ofS
¯
j
m
t is consistent with S
¯
11, and we can get the maximum value of S
¯
j
m
t by means of solving derivative of S
¯
j
m
t with respect to φ
m. Because the sine function term has nothing to do with φ
m, solving derivative of S
¯
j
m
t becomes solving derivative of S
¯
11.LetX
=
-
α
cos
⁡
θ
s
cos
⁡
φ
m
+
α
-
cos
⁡
θ
s
cos
⁡
φ
s, Y
=
-
α
sin
θ
s
cos
⁡
φ
m
+
α
-
sin
θ
s
cos
⁡
φ
s. Equation (30) becomes
(32)
S
¯
11
=
cos
⁡
X
+
cos
⁡
2
X
+
Y
2
+
cos
⁡
Y
+
cos
⁡
2
X
-
Y
2
.We have knownr
0
≪
R
j, r
0
≪
R
s, α
≪
1, α
-
≪
1, so X
≪
1, Y
≪
1. And then (32) can be given by
(33)
S
¯
11
=
1
-
X
2
2
+
1
-
2
X
+
Y
2
/
2
2
+
1
-
Y
2
2
+
1
-
2
X
-
Y
2
/
2
2
=
4
-
X
2
-
Y
2
.Solving the derivative ofS
¯
11 with respect to φ
m and setting d
S
¯
11
/
d
φ
m
=
0, we can get
(34)
cos
⁡
φ
m
=
α
-
cos
⁡
φ
s
α
.Because ofα
-
/
α
=
r
0
1
-
r
0
/
R
s
2
/
2
/
c
/
r
0
1
-
r
0
/
R
j
2
/
2
/
c
≈
1, we have
(35)
cos
⁡
φ
m
=
cos
⁡
φ
s
.From the above mathematical proof, we can know that whenφ
m
=
φ
s, S
¯
j
m
t 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 S
¯
j
m
t is only relevant with φ
m and has nothing to do with what R
j is. So the result obtained by scanning along with any circular arc whose radius is R
j
(
R
j
≫
r
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.Figure 5
Simplified search area of elevation angle.After the above analyses, the simplifying process about search area is finished, and the simplified search area only contains a circular line (see Figure3) 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.


## 3. Algorithm of VTR Direction Finding Method

### 3.1. Time Delay Correction

If we want to apply VTR direction finding method in practical engineering, the time delay correction is necessary. In practical engineering, the conversion rate of analog to digital converter (ADC) is limited; thus in order to improve the sampling accuracy, the radio frequency (RF) signal is usually converted to intermediate frequency (IF) signal via frequency conversion.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:
(36)
φ
=
2
π
f
c
τ
=
2
π
f
i
τ
′
.Thus,(37)
τ
′
=
f
c
τ
f
i
=
K
-
τ
,
where φ is the phase of input signal, f
c is the RF of input signal, f
i is the IF after frequency conversion, and K
-
=
f
c
/
f
i is the corrective factor of time delay.From (37) we can see that K
- is infinity if f
i
=
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.


### 3.2. VTR Algorithm for Estimating DOA

We have estimated the azimuth angle and elevation angle in Sections2.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 inT 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:
(38)
d
k
i
l
=
τ
P
k
i
l
f
s
,
where d
k
i
l denotes translated sample points in search area of azimuth angle and f
s is the sample rate.
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
(39)
d
k
i
l
=
τ
P
k
i
l
′
f
s
,
where τ
P
k
i
l
′ 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:(40)
d
k
j
m
=
τ
P
k
j
m
f
s
,
where d
k
j
m 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
(41)
d
k
j
m
=
τ
P
k
j
m
′
f
s
,
where τ
P
k
j
m
′ is the corrective time delay.
The time-reversed signal in Step2 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.


## 3.1. Time Delay Correction

If we want to apply VTR direction finding method in practical engineering, the time delay correction is necessary. In practical engineering, the conversion rate of analog to digital converter (ADC) is limited; thus in order to improve the sampling accuracy, the radio frequency (RF) signal is usually converted to intermediate frequency (IF) signal via frequency conversion.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:
(36)
φ
=
2
π
f
c
τ
=
2
π
f
i
τ
′
.Thus,(37)
τ
′
=
f
c
τ
f
i
=
K
-
τ
,
where φ is the phase of input signal, f
c is the RF of input signal, f
i is the IF after frequency conversion, and K
-
=
f
c
/
f
i is the corrective factor of time delay.From (37) we can see that K
- is infinity if f
i
=
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.


## 3.2. VTR Algorithm for Estimating DOA

We have estimated the azimuth angle and elevation angle in Sections2.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 inT 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:
(38)
d
k
i
l
=
τ
P
k
i
l
f
s
,
where d
k
i
l denotes translated sample points in search area of azimuth angle and f
s is the sample rate.
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
(39)
d
k
i
l
=
τ
P
k
i
l
′
f
s
,
where τ
P
k
i
l
′ 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:(40)
d
k
j
m
=
τ
P
k
j
m
f
s
,
where d
k
j
m 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
(41)
d
k
j
m
=
τ
P
k
j
m
′
f
s
,
where τ
P
k
j
m
′ is the corrective time delay.
The time-reversed signal in Step2 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.


## 4. Computer Simulations and Performance Analyses

Consider a UCA as shown in Figure1, 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 isS (20 km, 22°, 7.28 km), and the elevation angle is φ
s
=
2
0
∘
<
4
5
∘. 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
P
S
/
N
0
/
2
·
BW, P
S is the signal power, N
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 isf
i
=
700 kHz. The length of window for recording signal is T
=
4 ms and the centre of window is slidable.
(4)
Parameters of VTR direction finding process: set the processing bandwidth to be 40 kHz, and the corrective factor of time delay isK
-
=
10.7
×
1
0
6
/
700
×
1
0
3
=
10.7
/
0.7
=
15.286.
(5)
Parameters of DOA searching: the radius of circular line isr
i
=
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.Figure 6
Energy distribution for estimating azimuth angler
i
=
20 km.In order to observe the estimation error of azimuth angle in the whole search area, we set the radius of circular line to ber
i
=
i
Δ
r (Δ
r
=
500 m, i
=
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.Figure 7
RMSEs of azimuth angle for different radiuses.In the second experiment we study the estimation error of azimuth angle for different SNR levels and different record timeT. We consider the radius of circular line is r
i
=
20 km, the record time is T
=
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 T, we consider the radius of scanning circle as r
i
=
20 km 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 T increasing, and the estimation accuracy of VTR algorithm is better than compared methods at SNR = 0 dB.Figure 8
RMSEs of azimuth angle for different SNR levels.Figure 9
RMSEs of azimuth angle different record times.In the last experiment we study the estimation accuracy of elevation angle. Consider the parameters of source, channel, received signal, and VTR process which are the same as the first experiment. Set the radius of searching circular arc for elevation angle estimation to beR
j
=
21 km, and the searching increment is Δ
φ
=
1
∘. Figure 10 plots the RMSEs of elevation angle obtained by proposed method and compared methods at different SNR levels. We can see from Figure 10 that the RMSEs obtained by proposed method are less than compared methods at low SNR. And the estimation accuracy of proposed method is close to spectral MUSIC and RB-MUSIC when SNR > 12 dB. Figure 11 plots the RMSEs of elevation angle obtained by VTR algorithm and compared methods for different record time at SNR = 0 dB. It indicates that VTR algorithm delivers more accurate elevation angle estimation than compared methods for weak signal.Figure 10
RMSEs of elevation angle for different SNR levels.Figure 11
RMSEs of elevation angle for different record times.From Figures8, 9, 10, and 11 we can see that the VTR direction finding method can estimate 2D DOA accurately at low SNR levels. But the range of estimating elevation angle is restricted, and the proposed method is only suitable to detect the low altitude target whose elevation angle is φ
s
<
π
/
4.In order to evaluate the proposed method further, the computational complexity is discussed. As is known to all, the computational complexity of phase comparison method is very small, so here we only consider the computational complexities of spectral MUSIC, RB-MUSIC, and the proposed VTR method.For spectral MUSIC method, the computational complexity of algorithm mainly includes three parts: the covariance matrix calculation of received signals, the eigenvalue decomposition calculation, and the spectrum peak search. In spectral MUSIC, the computational complexity of the three parts isO
(
M
2
L
), O
(
M
3
), and O
(
M
2
+
(
M
2
+
M
)
(
N
1
×
N
2
)
). So the computational complexity of spectral MUSIC is O
(
M
2
L
+
M
3
+
M
2
+
(
M
2
+
M
)
(
N
1
×
N
2
)
), where, M is the number of elements in receive array, L is the number of samples, N
1 is the number of scanning points in azimuth angle, and N
2 is the number of scanning points in elevation angle. For RB-MUSIC, the computational complexity of algorithm mainly includes four parts: the covariance matrix calculation of received signals, the real beam space covariance matrix calculation, the eigenvalue decomposition calculation, and the spectrum peak search. In RB-MUSIC, the computational complexities of the given four parts are O
(
M
2
L
), O
(
B
M
2
+
B
2
M
), O
(
B
3
), and O
(
B
2
+
(
B
2
+
B
)
(
N
1
×
N
2
)
) separately, so the computational complexity of RB-MUSIC is O
(
M
2
L
+
B
M
2
+
B
2
M
+
B
3
+
B
2
+
(
B
2
+
B
)
(
N
1
×
N
2
)
), where B is the number of beam. For the proposed VTR method, the time delays for compensating signals can be calculated in advance if we choose the fixed r
i and R
j. Thus, the computational complexity of VTR method is only reflected in the energy calculation of scanning points for energy peak searching. If there only exists one processing unit for energy peak searching, the computational complexity of VTR method is O
(
L
(
N
1
+
N
2
)
). However, in our method, the parallel processing of energy peak searching can be performed by making use of more than one processing unit, so the computational complexity of VTR method can be simplified as the computational complexity of a processing unit. Hence when the number of processing units for energy peak searching is (
N
1
+
N
2
), the computational complexity of VTR method is simplified as O
(
L
(
N
1
+
N
2
)
/
(
N
1
+
N
2
)
)
=
O
(
L
).From the above analysis we could get to the conclusion that the computational complexity of VTR method only includes one term, which is closely related with the number of samplesL. And the proposed method can get a smaller computational complexity than spectral MUSIC and RB-MUSIC.In this paper, we are focusing on presenting a VTR method to estimate the DOA of single communication station and giving the mathematical proof and the simulation verification. We can get from the principle and process of 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.


## 5. 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φ
s
<
4
5
∘. 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.



---

*Source: 613692-2015-02-23.xml*