# Passive Radar Source Localization Based on Cuckoo Search NVTDOA

Authors: Gang Niu; Jie Gao; Tai-hang Du
Journal: International Journal of Antennas and Propagation
 (2020)

---

## Abstract

In order to meet the requirement of passive radar source localization in electronic warfare, the concept of the virtual time differences of arrival (VTDOA) is proposed by taking advantage of the characteristics of the same UAV in different positions at different times and the periodic rotation of radar pulse signal. The VTDOAs are the passive localization information defined as the time differences of the radar pulse transmission from the radar position to different virtual receivers. Firstly, a nonlinear VTDOA (NVTDOA) localization model is constructed. Moreover, sufficient conditions for accurately calculating the periodic integers in the model are analyzed, and the observability conditions of the localization model determined are deduced. Secondly, the convergence solution of the NVTDOA localization equation is obtained by Cuckoo search algorithm; thus, passive radar source localization is realized. Finally, the performance of the proposed method is verified by comparing with the existing methods.


---

## Body

## 1. Introduction

Passive localization is a fundamental problem that usually accomplishes localization task using the radiation signals of the target without radiating electromagnetic wave signals outwards, which has been widely used in many fields such as radar [1], sonar [2], sensor network [3], and wireless communication [4]. The passive localization method belongs to a two-step localization system. In the first step, some signal parameters in the target source transmission are extracted, such as received signal strength [5], time difference of arrival (TDOA) [6, 7], and angle of arrival (AOA) [8, 9]. In the second step, the target position is estimated by solving a set of nonlinear localization equations.For a long time, the passive localization method based on TDOA has been an important research topic [10]. For noncooperative radiation sources, it is particularly important to study TDOA localization technology for the single station because the time of signal emission is unknown. There are three main calculation methods of TDOA for noncooperative radiation sources. The first is to take the rising edge of the first pulse as the arrival time [11]. The pulse amplitude jump points from scratch are chosen as the reference point of the method, so it is easier and more accurate to determine the reference point. However, the signal amplitude of the first pulse is small, and the estimation error becomes larger when the pulse loss is more serious or the amplitude is smaller. The second is to take the average arrival time of all received pulses as the arrival time of radar signals [12]. The method has high stability but high computational complexity. If the received pulses overlap, the estimated performance will deteriorate. The third is to select the arrival time of the pulse with the highest amplitude in the main lobe pulse as the receiving time of radar signal [13]. When the main lobe beam is narrow, the arrival time will have higher stability. But if the main lobe beam is wider, there will be more maximum amplitude values in the pulse, which will degrade the estimation performance. In order to overcome the shortcomings of the above TDOA estimation, a norm approximation method was proposed by Wan [14]. The approximate mathematical expression of the main lobe beam pattern of the radar antenna is obtained at the receiving end, and the time point corresponding to the maximum value is used as the arrival time of radar signal to complete the estimation of VTDOA for radar signal.With the more and more extensive application of the unmanned aerial vehicle (UAV) in battlefield, there are some related research studies on the algorithm of using a UAV receiver as moving single station to locate a radar emitter. A passive location method of pulse information sources in the battlefield by using a single UAV was proposed by Nickens and Darko [15]. The method takes the rising time of pulse signals received by the UAV as the signal arrival time and builds localization equations based on UAV’s movement trajectory. However, when the pulse amplitude of the source is small or far away, the judgment error of the receiver on the rising moment of the pulse will be larger, and the localization accuracy will be worse. Kim et al. [16] proposed a way for a single drone to use a known database to provide rough location services for targets in the battlefield when the satellite positioning system is not working properly. To improve the performance of the location system in scenarios with a small number of receivers or when there are technical restrictions to deploy them in LoS condition, in reference [17–20], a localization algorithm based on time difference of arrival (TDOA) virtual sensors is proposed, where the reflection points are regarded as TDOA virtual sensors. And the algorithm is applied to localization of UAV in urban scenarios, ground vehicles, emitters, and RF emitters. Based on this, Hao et al. [21], Wan [14], and Zhu et al. [22] put forward a radar radiation source localization method based on passive synthetic aperture antenna array. The method utilizes the characteristics of a single-station small-load UAV in different positions at different times and the periodic rotation characteristics of radar pulse signals and uses VTDOA to complete the localization task. However, the observability of the localization model is not discussed, and the applicable conditions of the model are not analyzed in [14, 21, 22]. Moreover, it is assumed that the two radar signal waves received by the UAV are the same plane waves to obtain the closed solution of the positioning equation, and the radar rotation speed at the first received signal is used to calculate the radar antenna rotation time difference. Therefore, if the UAV runs at a high speed or the localization time difference is long, the localization accuracy of the target radar will become very poor.For the deficiency of the above algorithm, based on the motion characteristics of UAV and the periodic rotation characteristics of radar pulse signal, a new method of passive localization for the radar emitter based on Cuckoo search NVTDOA is proposed. In conclusion, the Cuckoo search NVTDOA method has made contributions in two aspects compared with the literature. First, the NVTDOA model is not sensitive to the speed of UAV or the localization time difference, so that it has more applicability. However, the approximate calculation of the rotation time difference in the VTDOA model makes its accuracy more sensitive to these, so VTDOA is only applicable to the positioning of UAV at low speed and short stroke. Second, the Cuckoo search algorithm is used instead of the closed-form solution algorithm to solve the localization equation, which finally improves the localization accuracy of the algorithm. In Section2, the traditional VTDOA localization model is introduced. In Section 3, first, the shortcomings of the VTDOA localization model are analyzed in the first paragraph, so an improved NVTDOA localization model is proposed considering the integral characteristics of radar antenna rotation time difference. Second, the NVTDOA localization model is described according to the traditional two-step localization method. The first step, the derivation of the NVTDOA localization model, is given in Section 3.1, and the second step, the estimation method of virtual time difference in the NVTDOA localization model, is introduced in Section 3.2. Therefore, it provides a localization model for the construction of the following localization equation. In Section 4, the observability of the NVTDOA model is introduced to verify the rationality of the model, and the theoretical basis is provided for the approximate calculation of the localization equationin the followi. Based on this, the construction and solution of the location equation are presented in Section 5. To realize the location of a radar emitter, the location models of two flight segments are established in Section 5.1, and the two models are combined to generate the location equations. According to the nonlinear characteristics of the equations, in Section 5.2, the Cuckoo search algorithm is proposed to solve the equations and finally realize the location of the radar emitter. In Section 6, VTDOA and NVTDOA are used to analyze the location of three different passive radars to verify the high accuracy and antinoise of the NVTDOA algorithm. In Section 7, the theoretical derivation and simulation are summarized and analyzed.


## 2. VTDOA Localization Model

The angular velocity of the target radar antenna rotating uniformly in counterclockwise direction is Ω, and the periodic rotation angle is Π. Taking theti moment of the UAV as the origin, a simple motion model of a UAV and the target radar is shown in Figure 1. The angle between the main lobe pulse direction of the far-field radiator and the UAV is φi. If the UAV moves in a straight line, it can be regarded as a constant.Figure 1
Simple motion model of a UAV and the target radar.Assuming that the radar rotation periodT=Π/Ω, which can be measured, the UAV receives the radar main lobe pulse signal at position Si at time ti, and after a period of time, the UAV moves to position Si+1 at time ti+1 and receives the radar main lobe pulse signal. Let P be the shortcut point of radiation source S0 from the UAV track. The length from S0 to P is D, and the length from Si to P is L. The distance from the initial point Si of the UAV to S0 is di and that from Si+1 to S0 is di+1. Similarly, the distance from Si+2 to S0 is di+2. The distance from location Si+1 to Si is li, that from Si+2 to Si+1 is li+1, and so on. Of these quantities, li and li+1 are known quantities. It can be seen from Figure 1 that, for the fixed radiation source S0, it can be characterized by X=D,LT. Then, the measurement equation of the variables to be solved can be written as follows:(1)Ψi+1,i=ti+1−ti=di+1−dic+Δψ+Ki+1,iT+Δti+1,i,where c=3.0×108m/s, which is the propagation speed of radar radio waves. Δti+1,i is the difference of TOA measurement noise between two locations of UAV, which is modeled as Gaussian white noise with zero mean value and standard deviation of σt. Ki+1,i represents the integer of radar rotation period, and it means that the time interval between two UAV locations is Ki+1,i radar period. Δψ is the rotation time difference of the radar antenna from Si to Si+1.In the VTDOA localization model, it is assumed that the linear velocity of the radar antenna during the rotation fromSi to Si+1 is the same as(2)vri=Ω⋅di.So(3)Δψ=lisinφivri=lisinφiΩdi.Then, the VTDOA localization model is that(4)VTDOAi+1,i=di+1−dic+lisinφiΩdi+Ki+1,iT+Δti+1,i.


## 3. NVTDOA Localization Model

### 3.1. NVTDOA Localization Model Based on Rotation Time Difference of Radar Antenna

In [14, 17, 18], in order to find the closed solution of the localization equation, it is assumed that the two radar signal waves received by UAV are the same plane waves. The rotation time difference of the radar antenna is calculated by using the radar rotation speed at the first received signal. Therefore, if the UAV runs at a high speed or the localization time difference is long, the localization accuracy of the target radar will become very poor. So, in this paper, an algorithm is presented to obtain the time difference of the radar antenna rotation based on definite integration of the rotation angle.As shown in Figure2, Δψ can be obtained by integrating the rotation angle θ of the radar antenna. The rotation time difference of the radar from Sθ to any small angle Δθ is Δψθ; then(5)Δψθ=ΔxΩdθ.Figure 2
The movement process of a UAV fromSi to Si+1.SoΔψ can be expressed as(6)Δψ=∫0θi+1,iΔψθ=∫0θi+1,iΔxΩdθ.According to the sine theorem,Δx can be obtained by(7)dθsinφ3=Δx/cosφ2sinΔθ,where φ2=π/2−φθ, φθ=φi+θ, φ3=π−Δθ−φ1, and φ1=φθ; then(8)Δx=dθcosφ2sinΔθsinφ3=dθcosπ/2−φi−θsinΔθsinπ−Δθ−φi−θ=dθsinφi+θsinΔθsinΔθ+φi+θ.Substituting equation (8) into equation (5) yields(9)Δψθ=dθsinφi+θsinΔθsinΔθ+φi+θΩdθ=sinφi+θsinΔθΩsinΔθ+φi+θ,where Δθ is an infinitesimal amount; then, sinΔθ≈Δθ, sinΔθ+φi+θ≈sinφi+θ, and equation (9) can be reduced to(10)Δψθ=ΔθΩ.Substituting equation (10) into equation (6), we have(11)Δψ=∫0θi+1,iΔθΩ=∫0θi+1,i1Ωdθ,where θi+1,i=cos−1di2+di+12−li2/2didi+1; then(12)Δψ=θi+1,iΩ=cos−1di2+di+12−li2/2didi+1Ω.Substituting equation (12) into equation (1) yields(13)Ψi+1,i=ti+1−ti=di+1−dic+θi+1,iΩ+Ki+1,iT+Δti+1,i.


### 3.2. Estimation Method of NVTDOA

When the radar system conducts periodic scanning, assuming that the UAV starts to collect the data at timet0, the UAV will not always be within the range of radar antenna fan beam irradiation. At the ith rotation period of radar, with the antenna rotating, the UAV will enter the range of the antenna main lobe fan beam at time ti and leave the main lobe irradiation area after receiving NP radar main lobe pulses continuously. And the center time of the kth main lobe pulse received in the ith rotation period is ti,k. Under the influence of beam gain directionality Gt, NP received pulse signals have different amplitudes. So the baseband signal waveform of the radar radiation source pulse signal received by UAV is shown in Figure 1, and PRI=1/PRF which is the pulse time interval, where PRF is the radar pulse repetition frequency. So the baseband signal of the linear frequency modulated (LFM) radar pulse received by the UAV can be expressed as follows [22]:(14)xt=PT⋅Gt⋅GR⋅at⋅ejφc+ξt,where ξt is additive random noise, PT is the transmitting power of the radar antenna, GR is the transceiver gains, and φc is the initial phase of the signal.As shown in Figure3, the aircraft will receive NP main lobe pulses in each rotation cycle of the radar. It is necessary to mine the hidden useful information from the received signal for radar location. The hidden information used in the paper is defined as NVTDOA. NVTDOA is defined as the difference of transmission time required for the radar signal to reach the same receiver at two different positions. In this paper, NVTDOA is represented by Ψj,i, and ti and tj represent the initial time for the same UAV to enter the ith and jth scanning periods of the main lobe of the radar antenna, respectively. So(15)Ψj,i=tj−ti.Figure 3
Received signal by UAV.In practical application, it is difficult to accurately estimateti and tj by xt of equation (14). As shown in Figure 3, let Ψj,i=tj,max−ti,max, where ti,max and tj,max are the peak times of Gt of the ith and jth rotation period, respectively. In order to improve the estimation accuracy of NVTDOA, the estimation method of ti,max and tj,max in this paper is referred to [21], which will not be described here.


## 3.1. NVTDOA Localization Model Based on Rotation Time Difference of Radar Antenna

In [14, 17, 18], in order to find the closed solution of the localization equation, it is assumed that the two radar signal waves received by UAV are the same plane waves. The rotation time difference of the radar antenna is calculated by using the radar rotation speed at the first received signal. Therefore, if the UAV runs at a high speed or the localization time difference is long, the localization accuracy of the target radar will become very poor. So, in this paper, an algorithm is presented to obtain the time difference of the radar antenna rotation based on definite integration of the rotation angle.As shown in Figure2, Δψ can be obtained by integrating the rotation angle θ of the radar antenna. The rotation time difference of the radar from Sθ to any small angle Δθ is Δψθ; then(5)Δψθ=ΔxΩdθ.Figure 2
The movement process of a UAV fromSi to Si+1.SoΔψ can be expressed as(6)Δψ=∫0θi+1,iΔψθ=∫0θi+1,iΔxΩdθ.According to the sine theorem,Δx can be obtained by(7)dθsinφ3=Δx/cosφ2sinΔθ,where φ2=π/2−φθ, φθ=φi+θ, φ3=π−Δθ−φ1, and φ1=φθ; then(8)Δx=dθcosφ2sinΔθsinφ3=dθcosπ/2−φi−θsinΔθsinπ−Δθ−φi−θ=dθsinφi+θsinΔθsinΔθ+φi+θ.Substituting equation (8) into equation (5) yields(9)Δψθ=dθsinφi+θsinΔθsinΔθ+φi+θΩdθ=sinφi+θsinΔθΩsinΔθ+φi+θ,where Δθ is an infinitesimal amount; then, sinΔθ≈Δθ, sinΔθ+φi+θ≈sinφi+θ, and equation (9) can be reduced to(10)Δψθ=ΔθΩ.Substituting equation (10) into equation (6), we have(11)Δψ=∫0θi+1,iΔθΩ=∫0θi+1,i1Ωdθ,where θi+1,i=cos−1di2+di+12−li2/2didi+1; then(12)Δψ=θi+1,iΩ=cos−1di2+di+12−li2/2didi+1Ω.Substituting equation (12) into equation (1) yields(13)Ψi+1,i=ti+1−ti=di+1−dic+θi+1,iΩ+Ki+1,iT+Δti+1,i.


## 3.2. Estimation Method of NVTDOA

When the radar system conducts periodic scanning, assuming that the UAV starts to collect the data at timet0, the UAV will not always be within the range of radar antenna fan beam irradiation. At the ith rotation period of radar, with the antenna rotating, the UAV will enter the range of the antenna main lobe fan beam at time ti and leave the main lobe irradiation area after receiving NP radar main lobe pulses continuously. And the center time of the kth main lobe pulse received in the ith rotation period is ti,k. Under the influence of beam gain directionality Gt, NP received pulse signals have different amplitudes. So the baseband signal waveform of the radar radiation source pulse signal received by UAV is shown in Figure 1, and PRI=1/PRF which is the pulse time interval, where PRF is the radar pulse repetition frequency. So the baseband signal of the linear frequency modulated (LFM) radar pulse received by the UAV can be expressed as follows [22]:(14)xt=PT⋅Gt⋅GR⋅at⋅ejφc+ξt,where ξt is additive random noise, PT is the transmitting power of the radar antenna, GR is the transceiver gains, and φc is the initial phase of the signal.As shown in Figure3, the aircraft will receive NP main lobe pulses in each rotation cycle of the radar. It is necessary to mine the hidden useful information from the received signal for radar location. The hidden information used in the paper is defined as NVTDOA. NVTDOA is defined as the difference of transmission time required for the radar signal to reach the same receiver at two different positions. In this paper, NVTDOA is represented by Ψj,i, and ti and tj represent the initial time for the same UAV to enter the ith and jth scanning periods of the main lobe of the radar antenna, respectively. So(15)Ψj,i=tj−ti.Figure 3
Received signal by UAV.In practical application, it is difficult to accurately estimateti and tj by xt of equation (14). As shown in Figure 3, let Ψj,i=tj,max−ti,max, where ti,max and tj,max are the peak times of Gt of the ith and jth rotation period, respectively. In order to improve the estimation accuracy of NVTDOA, the estimation method of ti,max and tj,max in this paper is referred to [21], which will not be described here.


## 4. Observability Analysis of Localization Model

The localization model is obtained in the previous section; the next problem is in which case the equation is solvable, namely, the observability problem of the model [23, 24]. In this paper, the observability conditions of the localization model are obtained by identifying the parameters of two systems.


### 4.1. Calculation Conditions for Integers of Radar Rotation Period

#### 4.1.1. Calculation Conditions with the Radar Rotation Period Known

When the radar rotation period is known, the period integer can be calculated by the following equation, takingSi and Si+1 as examples:(16)Ki+1,i=roundti+1−tiT,where round⋅ denotes the rounding integer operation. Substituting equation (13) into equation (16) yields(17)Ki+1,i=rounddi+1−dicT+θi+1,iΠ+Ki+1,i+Δti+1,iT.Sufficient conditions for accurate estimation which can be obtained by equation (17) are as follows:Condition (1a):di+1−di/cT+Δti+1,i/T≈0Condition (2a):θi+1,i<0.5ΠCondition (1a) has two meanings. The first is thatdi+1−di≪cT. The radar rotation period is generally in the order of seconds, so condition (1a) is evidently valid. The second is that Δti+1,i≪T, which is usually satisfied in the radar because the TOA measurement noise is typically on the order of ns and the radar radiation source period is typically on the order of s.The meaning of condition (2a) is that the angleθi+1,i between two measurements of the UAV and target is less than 0.5Π, so that the integer operation can be avoided to estimate Ki+1,i=±i, where i is an integer. The rotation angle period of most target radars is 2π; that is, Π=2π. When a UAV does not move radially relative to the target, we will have that θi+1,i<π=0.5Π, which satisfies the conditions.It can be seen from conditions (1a) and (2a) that when the radar rotation period is known,K can be accurately estimated.


#### 4.1.2. Calculation Conditions with the Radar Rotation Period Estimated

The radar antenna rotation period is unknown but can be estimated by the received radiation source pulse. It can be recorded asT^=T+ΔT, and ΔT is the estimation error. We have(18)Ki+1,i=roundti+1−tiT^.The integral fraction expansion is as follows:(19)ti+1−tiT^=di+1−dicT+ΔT+θi+1,iΩT+ΔT+Ki+1,iTT+ΔT+Δti+1,iT+ΔT=11+ΔT/Tdi+1−dicT+θi+1,iΠ+Ki+1,i+Δti+1,iT.According to equation (19), the sufficient conditions for accurately estimating Ki+1,i are as follows:Condition (1b):di+1−di/cT+Δti+1,i/T≈0Condition (2b):θi+1,i<0.5ΠCondition (3b):ΔT/T≈0The meaning of conditions (1b) and (2b) is consistent with conditions (1a) and (2a); the meaning of condition (3b) is that the estimation deviation of the rotation period of the radar antenna is negligible compared to the period itself, which is also the basic condition for the accuracy of the period estimation.


### 4.2. Observability Conditions of the Localization Model

The NVTDOA equation of the target radar position is nonlinear. In most papers, the observability theorem of Lee and Dum [25] is used to solve the observability condition of the equation. When the theorem is applied to discrete nonlinear systems, the equivalent conclusion is to judge the rank of the Gramer matrix. For a nonlinear system as shown in equation (20), let the system set be S and the Gramer matrix of the state variables xi0∗ in the set is as in equation (21):(20)xn=fxn,n;x1=x0,yn=hxn,n,(21)Γi0,i0+N−1=Hi0,Hi0+1Φi0,i0+1,…,Hi0+N−1Φi0,i0+N−1T,where Hi=∂yn/∂xn, which is the Jacobian matrix, and Φ is a state transition matrix. If there is an integer N such that rankΓi,i+N−1=n (state dimension), then the system is completely observable for S [26].


#### 4.2.1. Observability Analysis with the Radar Rotation Period Known

When the radar period is known, the state transition matrixΦi,i+1=I2×2 and the measurement equation of the variable X=D,LT can be simplified as(22)NVTDOAi+1,i=di+1−dic+θi+1,iΩ+Ki+1,iT+Δti+1,i.TheJacobi matrix of NVTDOAj,i at X is(23)Hi=∂NVTDOAi+1,i∂X=∂NVTDOAi+1,i∂D,∂NVTDOAi+1,i∂LΤ.From Figure1, it can be seen that di=D2+L2anddi+1=D2+L−li2. By substituting θi+1,i=tan−1L/D−tan−1L−li/D into equation (23), we can get(24)Hi=∂NVTDOAi+1,i∂X=1cDdi+1−Ddi+1ΩL−lidi+12−Ldi21cL−lidi+1−Ldi+1ΩDdi2−Ddi+12Τ.Under far-field radiation source settings, the radar source is a two-dimensional target. Let the position of the UAV receiving the radiation source pulse beSi, Si+1, and Si+2, and Gramer matrix based on NVTDOA by three measurements is(25)Γi,i+2=Hi,Hi+1Τ,where(26)Hi+1=∂NVTDOAi+2,i+1∂X=1cDdi+2−Ddi+1+1ΩL−li−li+1di+22−L−lidi+121cL−li−li+1di+2−L−lidi+1+1ΩDdi+12−Ddi+22Τ.The determinant of equation (25) is that(27)detΓi,i+2=Dc2M1didi+1di+2+D2cΩM2di2di+12di+22+1cΩM3di2di+12di+22+DΩ2M4di2di+12di+22,where(28)M1=li+1di+1−di−lidi+2−di+1,M2=di2di+2−di+1−di+12di+2−di+di+22di+1−di,M3=di2L−li−li+1L−lidi+2−di+1−di+12LL−li−li+1di+2−di+di+22LL−lidi+1−di,M4=di+12−di+22li+di+12−di2li+1.Discussions aboutMi are as follows:(1)
Ifli=0, then M1≠0 when li+1≠0 and di≠di+1. But li=0 and di≠di+1 are contradictory, so this condition is not valid(2)
Ifli+1=0, then M1≠0 when li≠0 and di+1≠di+2. But li+1=0 and di+1≠di+2 are contradictory, so this condition is not validTherefore, from (1) and (2), we can know M1≠0 when li≠0, li+1≠0, and di≠di+1≠di+2(3)
Whendi≠di+1≠di+2≠0, M2≠0(4)
Ifli=0, then di=di+1. So M3=0. If li+1=0, then di+1=di+2. So M3=0. Therefore, when li≠0,li+1≠0, and di≠di+1≠di+2, M3≠0(5)
Whendi≠di+1≠di+2≠0, li≠0, and li+1≠0, M4≠0From (1)–(5) and (27), it can be seen that when D≠0, di≠di+1≠di+2≠0, li≠0, and li+1≠0, then detΓi,i+2≠0. It shows that D and L can be estimated by three NVTDOA measurements when UAV must move but not radially. So the target radar position can be determined.


#### 4.2.2. Observability Analysis with the Radar Rotation Period Estimated

When the periodT is unknown, the period estimate X∗=D,L,TT and the state transition matrix Φi,i+1=I3×3 need to be added to analyze the observability. Without affecting observability, the period integer of two adjacent measurements of UAV is assumed to be 1:(29)Hi∗=∂NVTDOAi+3,i+2∂X∗=1cDdi+1−Ddi+1ΩL−lidi+12−Ldi21cL−lidi+1−Ldi+1ΩDdi2−Ddi+121Τ,Hi+1∗=∂NVTDOAi+3,i+2∂X∗=1cDdi+2−Ddi+1+1ΩL−−li−li+1di+22−L−lidi+121cL−li−li+1di+2−L−lidi+1+1ΩDdi+12−Ddi+221Τ,Hi+2∗=∂NVTDOAi+3,i+2∂X∗=1cDdi+3−Ddi+2+1ΩL−li−li+1−li+2di+32−L−li−li+1di+221cL−li−li+1−li+2di+3−L−li−li+1di+2+1ΩDdi+22−Ddi+321Τ.Similarly, the Gramer matrix based on NVTDOA by four measurements is(30)Γ∗i,i+3=Hi∗,Hi+1∗,Hi+2∗Τ,(31)detΓ∗i,i+3=detHi+1,Hi+2−detHi,Hi+2+detHi,Hi+1.From the calculation in the previous section, we can see the following:(1)
WhenD≠0, di≠di+1≠di+2≠0, li≠0, and li+1≠0, then detHi,Hi+1≠0(2)
WhenD≠0, di+1≠di+2≠di+3≠0, li+1≠0, and li+2≠0, then detHi+1,Hi+2≠0(3)
WhenD≠0, di≠di+2≠di+3≠0, li+1+li≠0, and li+2≠0, then detHi,Hi+2≠0In summary, the observability conditions of equation (30) can be obtained; that is, D≠0 and di≠di+1≠di+2≠di+3≠0, li≠0,li+1≠0,li+2≠0. That is to say, D, L, and T can be estimated by four NVTDOA measurements when UAV must move but not radially. So the target radar position can be determined.


## 4.1. Calculation Conditions for Integers of Radar Rotation Period

### 4.1.1. Calculation Conditions with the Radar Rotation Period Known

When the radar rotation period is known, the period integer can be calculated by the following equation, takingSi and Si+1 as examples:(16)Ki+1,i=roundti+1−tiT,where round⋅ denotes the rounding integer operation. Substituting equation (13) into equation (16) yields(17)Ki+1,i=rounddi+1−dicT+θi+1,iΠ+Ki+1,i+Δti+1,iT.Sufficient conditions for accurate estimation which can be obtained by equation (17) are as follows:Condition (1a):di+1−di/cT+Δti+1,i/T≈0Condition (2a):θi+1,i<0.5ΠCondition (1a) has two meanings. The first is thatdi+1−di≪cT. The radar rotation period is generally in the order of seconds, so condition (1a) is evidently valid. The second is that Δti+1,i≪T, which is usually satisfied in the radar because the TOA measurement noise is typically on the order of ns and the radar radiation source period is typically on the order of s.The meaning of condition (2a) is that the angleθi+1,i between two measurements of the UAV and target is less than 0.5Π, so that the integer operation can be avoided to estimate Ki+1,i=±i, where i is an integer. The rotation angle period of most target radars is 2π; that is, Π=2π. When a UAV does not move radially relative to the target, we will have that θi+1,i<π=0.5Π, which satisfies the conditions.It can be seen from conditions (1a) and (2a) that when the radar rotation period is known,K can be accurately estimated.


### 4.1.2. Calculation Conditions with the Radar Rotation Period Estimated

The radar antenna rotation period is unknown but can be estimated by the received radiation source pulse. It can be recorded asT^=T+ΔT, and ΔT is the estimation error. We have(18)Ki+1,i=roundti+1−tiT^.The integral fraction expansion is as follows:(19)ti+1−tiT^=di+1−dicT+ΔT+θi+1,iΩT+ΔT+Ki+1,iTT+ΔT+Δti+1,iT+ΔT=11+ΔT/Tdi+1−dicT+θi+1,iΠ+Ki+1,i+Δti+1,iT.According to equation (19), the sufficient conditions for accurately estimating Ki+1,i are as follows:Condition (1b):di+1−di/cT+Δti+1,i/T≈0Condition (2b):θi+1,i<0.5ΠCondition (3b):ΔT/T≈0The meaning of conditions (1b) and (2b) is consistent with conditions (1a) and (2a); the meaning of condition (3b) is that the estimation deviation of the rotation period of the radar antenna is negligible compared to the period itself, which is also the basic condition for the accuracy of the period estimation.


## 4.1.1. Calculation Conditions with the Radar Rotation Period Known

When the radar rotation period is known, the period integer can be calculated by the following equation, takingSi and Si+1 as examples:(16)Ki+1,i=roundti+1−tiT,where round⋅ denotes the rounding integer operation. Substituting equation (13) into equation (16) yields(17)Ki+1,i=rounddi+1−dicT+θi+1,iΠ+Ki+1,i+Δti+1,iT.Sufficient conditions for accurate estimation which can be obtained by equation (17) are as follows:Condition (1a):di+1−di/cT+Δti+1,i/T≈0Condition (2a):θi+1,i<0.5ΠCondition (1a) has two meanings. The first is thatdi+1−di≪cT. The radar rotation period is generally in the order of seconds, so condition (1a) is evidently valid. The second is that Δti+1,i≪T, which is usually satisfied in the radar because the TOA measurement noise is typically on the order of ns and the radar radiation source period is typically on the order of s.The meaning of condition (2a) is that the angleθi+1,i between two measurements of the UAV and target is less than 0.5Π, so that the integer operation can be avoided to estimate Ki+1,i=±i, where i is an integer. The rotation angle period of most target radars is 2π; that is, Π=2π. When a UAV does not move radially relative to the target, we will have that θi+1,i<π=0.5Π, which satisfies the conditions.It can be seen from conditions (1a) and (2a) that when the radar rotation period is known,K can be accurately estimated.


## 4.1.2. Calculation Conditions with the Radar Rotation Period Estimated

The radar antenna rotation period is unknown but can be estimated by the received radiation source pulse. It can be recorded asT^=T+ΔT, and ΔT is the estimation error. We have(18)Ki+1,i=roundti+1−tiT^.The integral fraction expansion is as follows:(19)ti+1−tiT^=di+1−dicT+ΔT+θi+1,iΩT+ΔT+Ki+1,iTT+ΔT+Δti+1,iT+ΔT=11+ΔT/Tdi+1−dicT+θi+1,iΠ+Ki+1,i+Δti+1,iT.According to equation (19), the sufficient conditions for accurately estimating Ki+1,i are as follows:Condition (1b):di+1−di/cT+Δti+1,i/T≈0Condition (2b):θi+1,i<0.5ΠCondition (3b):ΔT/T≈0The meaning of conditions (1b) and (2b) is consistent with conditions (1a) and (2a); the meaning of condition (3b) is that the estimation deviation of the rotation period of the radar antenna is negligible compared to the period itself, which is also the basic condition for the accuracy of the period estimation.


## 4.2. Observability Conditions of the Localization Model

The NVTDOA equation of the target radar position is nonlinear. In most papers, the observability theorem of Lee and Dum [25] is used to solve the observability condition of the equation. When the theorem is applied to discrete nonlinear systems, the equivalent conclusion is to judge the rank of the Gramer matrix. For a nonlinear system as shown in equation (20), let the system set be S and the Gramer matrix of the state variables xi0∗ in the set is as in equation (21):(20)xn=fxn,n;x1=x0,yn=hxn,n,(21)Γi0,i0+N−1=Hi0,Hi0+1Φi0,i0+1,…,Hi0+N−1Φi0,i0+N−1T,where Hi=∂yn/∂xn, which is the Jacobian matrix, and Φ is a state transition matrix. If there is an integer N such that rankΓi,i+N−1=n (state dimension), then the system is completely observable for S [26].


### 4.2.1. Observability Analysis with the Radar Rotation Period Known

When the radar period is known, the state transition matrixΦi,i+1=I2×2 and the measurement equation of the variable X=D,LT can be simplified as(22)NVTDOAi+1,i=di+1−dic+θi+1,iΩ+Ki+1,iT+Δti+1,i.TheJacobi matrix of NVTDOAj,i at X is(23)Hi=∂NVTDOAi+1,i∂X=∂NVTDOAi+1,i∂D,∂NVTDOAi+1,i∂LΤ.From Figure1, it can be seen that di=D2+L2anddi+1=D2+L−li2. By substituting θi+1,i=tan−1L/D−tan−1L−li/D into equation (23), we can get(24)Hi=∂NVTDOAi+1,i∂X=1cDdi+1−Ddi+1ΩL−lidi+12−Ldi21cL−lidi+1−Ldi+1ΩDdi2−Ddi+12Τ.Under far-field radiation source settings, the radar source is a two-dimensional target. Let the position of the UAV receiving the radiation source pulse beSi, Si+1, and Si+2, and Gramer matrix based on NVTDOA by three measurements is(25)Γi,i+2=Hi,Hi+1Τ,where(26)Hi+1=∂NVTDOAi+2,i+1∂X=1cDdi+2−Ddi+1+1ΩL−li−li+1di+22−L−lidi+121cL−li−li+1di+2−L−lidi+1+1ΩDdi+12−Ddi+22Τ.The determinant of equation (25) is that(27)detΓi,i+2=Dc2M1didi+1di+2+D2cΩM2di2di+12di+22+1cΩM3di2di+12di+22+DΩ2M4di2di+12di+22,where(28)M1=li+1di+1−di−lidi+2−di+1,M2=di2di+2−di+1−di+12di+2−di+di+22di+1−di,M3=di2L−li−li+1L−lidi+2−di+1−di+12LL−li−li+1di+2−di+di+22LL−lidi+1−di,M4=di+12−di+22li+di+12−di2li+1.Discussions aboutMi are as follows:(1)
Ifli=0, then M1≠0 when li+1≠0 and di≠di+1. But li=0 and di≠di+1 are contradictory, so this condition is not valid(2)
Ifli+1=0, then M1≠0 when li≠0 and di+1≠di+2. But li+1=0 and di+1≠di+2 are contradictory, so this condition is not validTherefore, from (1) and (2), we can know M1≠0 when li≠0, li+1≠0, and di≠di+1≠di+2(3)
Whendi≠di+1≠di+2≠0, M2≠0(4)
Ifli=0, then di=di+1. So M3=0. If li+1=0, then di+1=di+2. So M3=0. Therefore, when li≠0,li+1≠0, and di≠di+1≠di+2, M3≠0(5)
Whendi≠di+1≠di+2≠0, li≠0, and li+1≠0, M4≠0From (1)–(5) and (27), it can be seen that when D≠0, di≠di+1≠di+2≠0, li≠0, and li+1≠0, then detΓi,i+2≠0. It shows that D and L can be estimated by three NVTDOA measurements when UAV must move but not radially. So the target radar position can be determined.


### 4.2.2. Observability Analysis with the Radar Rotation Period Estimated

When the periodT is unknown, the period estimate X∗=D,L,TT and the state transition matrix Φi,i+1=I3×3 need to be added to analyze the observability. Without affecting observability, the period integer of two adjacent measurements of UAV is assumed to be 1:(29)Hi∗=∂NVTDOAi+3,i+2∂X∗=1cDdi+1−Ddi+1ΩL−lidi+12−Ldi21cL−lidi+1−Ldi+1ΩDdi2−Ddi+121Τ,Hi+1∗=∂NVTDOAi+3,i+2∂X∗=1cDdi+2−Ddi+1+1ΩL−−li−li+1di+22−L−lidi+121cL−li−li+1di+2−L−lidi+1+1ΩDdi+12−Ddi+221Τ,Hi+2∗=∂NVTDOAi+3,i+2∂X∗=1cDdi+3−Ddi+2+1ΩL−li−li+1−li+2di+32−L−li−li+1di+221cL−li−li+1−li+2di+3−L−li−li+1di+2+1ΩDdi+22−Ddi+321Τ.Similarly, the Gramer matrix based on NVTDOA by four measurements is(30)Γ∗i,i+3=Hi∗,Hi+1∗,Hi+2∗Τ,(31)detΓ∗i,i+3=detHi+1,Hi+2−detHi,Hi+2+detHi,Hi+1.From the calculation in the previous section, we can see the following:(1)
WhenD≠0, di≠di+1≠di+2≠0, li≠0, and li+1≠0, then detHi,Hi+1≠0(2)
WhenD≠0, di+1≠di+2≠di+3≠0, li+1≠0, and li+2≠0, then detHi+1,Hi+2≠0(3)
WhenD≠0, di≠di+2≠di+3≠0, li+1+li≠0, and li+2≠0, then detHi,Hi+2≠0In summary, the observability conditions of equation (30) can be obtained; that is, D≠0 and di≠di+1≠di+2≠di+3≠0, li≠0,li+1≠0,li+2≠0. That is to say, D, L, and T can be estimated by four NVTDOA measurements when UAV must move but not radially. So the target radar position can be determined.


## 4.2.1. Observability Analysis with the Radar Rotation Period Known

When the radar period is known, the state transition matrixΦi,i+1=I2×2 and the measurement equation of the variable X=D,LT can be simplified as(22)NVTDOAi+1,i=di+1−dic+θi+1,iΩ+Ki+1,iT+Δti+1,i.TheJacobi matrix of NVTDOAj,i at X is(23)Hi=∂NVTDOAi+1,i∂X=∂NVTDOAi+1,i∂D,∂NVTDOAi+1,i∂LΤ.From Figure1, it can be seen that di=D2+L2anddi+1=D2+L−li2. By substituting θi+1,i=tan−1L/D−tan−1L−li/D into equation (23), we can get(24)Hi=∂NVTDOAi+1,i∂X=1cDdi+1−Ddi+1ΩL−lidi+12−Ldi21cL−lidi+1−Ldi+1ΩDdi2−Ddi+12Τ.Under far-field radiation source settings, the radar source is a two-dimensional target. Let the position of the UAV receiving the radiation source pulse beSi, Si+1, and Si+2, and Gramer matrix based on NVTDOA by three measurements is(25)Γi,i+2=Hi,Hi+1Τ,where(26)Hi+1=∂NVTDOAi+2,i+1∂X=1cDdi+2−Ddi+1+1ΩL−li−li+1di+22−L−lidi+121cL−li−li+1di+2−L−lidi+1+1ΩDdi+12−Ddi+22Τ.The determinant of equation (25) is that(27)detΓi,i+2=Dc2M1didi+1di+2+D2cΩM2di2di+12di+22+1cΩM3di2di+12di+22+DΩ2M4di2di+12di+22,where(28)M1=li+1di+1−di−lidi+2−di+1,M2=di2di+2−di+1−di+12di+2−di+di+22di+1−di,M3=di2L−li−li+1L−lidi+2−di+1−di+12LL−li−li+1di+2−di+di+22LL−lidi+1−di,M4=di+12−di+22li+di+12−di2li+1.Discussions aboutMi are as follows:(1)
Ifli=0, then M1≠0 when li+1≠0 and di≠di+1. But li=0 and di≠di+1 are contradictory, so this condition is not valid(2)
Ifli+1=0, then M1≠0 when li≠0 and di+1≠di+2. But li+1=0 and di+1≠di+2 are contradictory, so this condition is not validTherefore, from (1) and (2), we can know M1≠0 when li≠0, li+1≠0, and di≠di+1≠di+2(3)
Whendi≠di+1≠di+2≠0, M2≠0(4)
Ifli=0, then di=di+1. So M3=0. If li+1=0, then di+1=di+2. So M3=0. Therefore, when li≠0,li+1≠0, and di≠di+1≠di+2, M3≠0(5)
Whendi≠di+1≠di+2≠0, li≠0, and li+1≠0, M4≠0From (1)–(5) and (27), it can be seen that when D≠0, di≠di+1≠di+2≠0, li≠0, and li+1≠0, then detΓi,i+2≠0. It shows that D and L can be estimated by three NVTDOA measurements when UAV must move but not radially. So the target radar position can be determined.


## 4.2.2. Observability Analysis with the Radar Rotation Period Estimated

When the periodT is unknown, the period estimate X∗=D,L,TT and the state transition matrix Φi,i+1=I3×3 need to be added to analyze the observability. Without affecting observability, the period integer of two adjacent measurements of UAV is assumed to be 1:(29)Hi∗=∂NVTDOAi+3,i+2∂X∗=1cDdi+1−Ddi+1ΩL−lidi+12−Ldi21cL−lidi+1−Ldi+1ΩDdi2−Ddi+121Τ,Hi+1∗=∂NVTDOAi+3,i+2∂X∗=1cDdi+2−Ddi+1+1ΩL−−li−li+1di+22−L−lidi+121cL−li−li+1di+2−L−lidi+1+1ΩDdi+12−Ddi+221Τ,Hi+2∗=∂NVTDOAi+3,i+2∂X∗=1cDdi+3−Ddi+2+1ΩL−li−li+1−li+2di+32−L−li−li+1di+221cL−li−li+1−li+2di+3−L−li−li+1di+2+1ΩDdi+22−Ddi+321Τ.Similarly, the Gramer matrix based on NVTDOA by four measurements is(30)Γ∗i,i+3=Hi∗,Hi+1∗,Hi+2∗Τ,(31)detΓ∗i,i+3=detHi+1,Hi+2−detHi,Hi+2+detHi,Hi+1.From the calculation in the previous section, we can see the following:(1)
WhenD≠0, di≠di+1≠di+2≠0, li≠0, and li+1≠0, then detHi,Hi+1≠0(2)
WhenD≠0, di+1≠di+2≠di+3≠0, li+1≠0, and li+2≠0, then detHi+1,Hi+2≠0(3)
WhenD≠0, di≠di+2≠di+3≠0, li+1+li≠0, and li+2≠0, then detHi,Hi+2≠0In summary, the observability conditions of equation (30) can be obtained; that is, D≠0 and di≠di+1≠di+2≠di+3≠0, li≠0,li+1≠0,li+2≠0. That is to say, D, L, and T can be estimated by four NVTDOA measurements when UAV must move but not radially. So the target radar position can be determined.


## 5. Construction and Solution of NVTDOA Localization Equation Based on Cuckoo Search

### 5.1. Construction of NVTDOA Localization Equation

When the UAV moves, we have(32)Ψj,i=tj−ti.The value ofΨj,i can be estimated by Section 3.2. Referring to equation (13), it can be concluded that(33)Ψj,i≈dj−dic+θj,iΩ+Kj,iT,where Kj,i=roundΨj,i/T. The distance from Si to S0 is di at ti time, the distance from Sj to S0 is dj, and the distance from Sj to Si is li.It should be noted that from the conclusion of Section4.2, at low speed the observability condition of the model is related to the change of the position of the UAV and has nothing to do with its speed. If the position of the UAV is recorded during the movement, the error will be caused because the initial time ti is not well determined. In order to locate Si accurately, the UAV hovers approximately before ti to obtain the exact position of the main lobe of the receiving radar antenna during the ith scanning period. After receiving the main lobe pulse in this period, the UAV flies to position Sj until receiving the main lobe pulse in the jth scanning period of the radar antenna.As shown in Figure1, the geometric relationship diagram of periodic emitter location based on NVTDOA is obtained. For convenience of representation, let Si+1=Sj and Si+2=Sl. In the Cartesian rectangular two-dimensional coordinate system, the radar emitter is located at S0=x0,y0T, UAV is located at Si=xi,yiT, Sj=xj,yjT, and Sl=xl,ylT. It is known from the Beidou localization system that its error can be ignored. Based on the accurate estimation of Ψj,i, equation (33) can be reduced to(34)Ψj,i≈dj−dic+θj,iΩ+Kj,iT≈θj,iΩ+Kj,iT,where θj,i=cos−1di2+dj2−li2/2didj, di=xi−x02+yi−y02, dj=xj−x02+yj−y02, li=xj−xi2+yj−yi2, and Ω=Π/T. The radar system rotation period angle is Π=360∘. The approximation method used in the equation is the same as condition (1a) in Section 4. Substituting di, dj, and li into θj,i, we have(35)θj,i=cos−1x0−xix0−xj+y0−yiy0−yjx0−xi2+y0−yi2⋅x0−xj2+y0−yj2.Substituting equation (35) into equation (34) yields(36)Ψj,i⋅Ω−Kj,i⋅Π−cos−1x0−xix0−xj+y0−yiy0−yjx0−xi2+y0−yi2⋅x0−xj2+y0−yj2=0.Equation (36) constructs a nonlinear equation based on the position S0=x0,y0T of the radar radiation source. Similarly, when the UAV flies from Si=xj,yjT to Sl=xl,ylT, there will be(37)Ψl,j⋅Ω−Kl,j⋅Π−cos−1x0−xjx0−xl+y0−yjy0−ylx0−xj2+y0−yj2⋅x0−xl2+y0−yl2=0.By solving the nonlinear equations consisting of equations (36) and (37), the localization results can be obtained that S^0=x^0,y^0T.


### 5.2. Solution of Nonlinear Localization Equations Based on Cuckoo Search Algorithm

#### 5.2.1. Cuckoo Search Algorithm

The cuckoo search (CS) algorithm is generated by the simulation study of cuckoo’s nesting behavior by Yang Xin-she and DEB Suash. The individuals in birds flocking select nests by Levy flight and abandon by discovered probabilityPa, so the birds population is updated. The paths and locations of cuckoo nests are updated as follows:(38)xit+1=xit+α0xit−xbestt⊕Levyλ,i=1,2,…,n,where xit=xi,1t,…,xi,dtT is the ith solution of the tth generation, α0=0.01, xbestt is the optimal solution of the tth generation, Levyλ∼ϕ×u/v1/λ is the dot product, n is the number of nests, and Levyλ∼ϕ×u/v1/λ is the Levy random search path whose moving step size obeys levy probability distribution as(39)Levyλ∼ϕ×uv1/λ,where u and v follow normal distribution and ϕ=Γ1+λsinλπ/2/Γ1+λ/2×21+λ/2λ, where Γ is the standard Gamma function.Combining equations (38) and (39), in CS algorithm, the new solution xit+1 is generated by(40)xit+1=xit+0.01xit−xbestt⊕ϕ×uv1/λ,i=1,2,…,n.After discarding a part of solutions according to the discovery probability, a new solution is generated by(41)xit+1=xit+rxjt−xet,where r is the random number of uniform distribution in (0, 1) interval and xjt and xet are random solutions of the tth generation.


#### 5.2.2. Solution of Nonlinear Localization Equations Based on Cuckoo Search Algorithm

The localization equations consisting of equations (36) and (37) are recorded as(42)f1x,y=0,f2x,y=0.When equation (42) is transformed into the unconstrained optimization problem, the objective function of the cuckoo search algorithm is(43)Fw=∑i=12fiw2,where w=x,yT. So solving equation (42) is to find the minimum solution of equation (43).To sum up, the solution process of CS is as follows:(1)
Initialization. The population number of the cuckoo algorithm is m, the maximum number of iterations is K, and the discovered probability is Pa. Moreover, m nest positions are randomly generated and the optimal nest position of the current generation,w0,bestbest∈1,2,…,m, and the optimal fitness value, Fmin, are found. In the simulation experiment of the paper, m = 20, K = 2000, and Pa = 0.25.(2)
Update of Levy Flight Path (Loop Body). The optimal value of the previous generation, wk−1,best, is reserved. K represents the current evolutionary generation. And a new nest location is generated by equation (40). Then we compare the optimal fitness value of the new nest location with that of the previous nest location and keep the solution with the smaller optimal fitness value which is expressed as Ik=wk,1,wk,2,…,wk,m to the next nest location.(3)
Update of the Discovered Probability (Loop Body). For each individual in qk=wk,1,wk,2,⋯,wk,m, r is randomly generated as the probability that the bird’s egg is found by the nest owner, that is, the discovered probability. If r < Pa, the nest value is unchanged; otherwise, the nest is updated according to equation (41). The new nest is compared with qk=wk,1,wk,2,⋯,wk,m, and then the better solution which is expressed as qk=wk,1,wk,2,…,wk,m is kept to the next generation.(4)
Judgment Body. The optimal nest position wk,best in qk and the fitness value Fmin of the solution are calculated. At the same time, it is determined whether the termination condition of the algorithm (the generation is equal to K) is satisfied. If it is satisfied, the global optimal solution is obtained. Otherwise, step 2 is returned, and the iterative operation continues until the termination condition is satisfied.


## 5.1. Construction of NVTDOA Localization Equation

When the UAV moves, we have(32)Ψj,i=tj−ti.The value ofΨj,i can be estimated by Section 3.2. Referring to equation (13), it can be concluded that(33)Ψj,i≈dj−dic+θj,iΩ+Kj,iT,where Kj,i=roundΨj,i/T. The distance from Si to S0 is di at ti time, the distance from Sj to S0 is dj, and the distance from Sj to Si is li.It should be noted that from the conclusion of Section4.2, at low speed the observability condition of the model is related to the change of the position of the UAV and has nothing to do with its speed. If the position of the UAV is recorded during the movement, the error will be caused because the initial time ti is not well determined. In order to locate Si accurately, the UAV hovers approximately before ti to obtain the exact position of the main lobe of the receiving radar antenna during the ith scanning period. After receiving the main lobe pulse in this period, the UAV flies to position Sj until receiving the main lobe pulse in the jth scanning period of the radar antenna.As shown in Figure1, the geometric relationship diagram of periodic emitter location based on NVTDOA is obtained. For convenience of representation, let Si+1=Sj and Si+2=Sl. In the Cartesian rectangular two-dimensional coordinate system, the radar emitter is located at S0=x0,y0T, UAV is located at Si=xi,yiT, Sj=xj,yjT, and Sl=xl,ylT. It is known from the Beidou localization system that its error can be ignored. Based on the accurate estimation of Ψj,i, equation (33) can be reduced to(34)Ψj,i≈dj−dic+θj,iΩ+Kj,iT≈θj,iΩ+Kj,iT,where θj,i=cos−1di2+dj2−li2/2didj, di=xi−x02+yi−y02, dj=xj−x02+yj−y02, li=xj−xi2+yj−yi2, and Ω=Π/T. The radar system rotation period angle is Π=360∘. The approximation method used in the equation is the same as condition (1a) in Section 4. Substituting di, dj, and li into θj,i, we have(35)θj,i=cos−1x0−xix0−xj+y0−yiy0−yjx0−xi2+y0−yi2⋅x0−xj2+y0−yj2.Substituting equation (35) into equation (34) yields(36)Ψj,i⋅Ω−Kj,i⋅Π−cos−1x0−xix0−xj+y0−yiy0−yjx0−xi2+y0−yi2⋅x0−xj2+y0−yj2=0.Equation (36) constructs a nonlinear equation based on the position S0=x0,y0T of the radar radiation source. Similarly, when the UAV flies from Si=xj,yjT to Sl=xl,ylT, there will be(37)Ψl,j⋅Ω−Kl,j⋅Π−cos−1x0−xjx0−xl+y0−yjy0−ylx0−xj2+y0−yj2⋅x0−xl2+y0−yl2=0.By solving the nonlinear equations consisting of equations (36) and (37), the localization results can be obtained that S^0=x^0,y^0T.


## 5.2. Solution of Nonlinear Localization Equations Based on Cuckoo Search Algorithm

### 5.2.1. Cuckoo Search Algorithm

The cuckoo search (CS) algorithm is generated by the simulation study of cuckoo’s nesting behavior by Yang Xin-she and DEB Suash. The individuals in birds flocking select nests by Levy flight and abandon by discovered probabilityPa, so the birds population is updated. The paths and locations of cuckoo nests are updated as follows:(38)xit+1=xit+α0xit−xbestt⊕Levyλ,i=1,2,…,n,where xit=xi,1t,…,xi,dtT is the ith solution of the tth generation, α0=0.01, xbestt is the optimal solution of the tth generation, Levyλ∼ϕ×u/v1/λ is the dot product, n is the number of nests, and Levyλ∼ϕ×u/v1/λ is the Levy random search path whose moving step size obeys levy probability distribution as(39)Levyλ∼ϕ×uv1/λ,where u and v follow normal distribution and ϕ=Γ1+λsinλπ/2/Γ1+λ/2×21+λ/2λ, where Γ is the standard Gamma function.Combining equations (38) and (39), in CS algorithm, the new solution xit+1 is generated by(40)xit+1=xit+0.01xit−xbestt⊕ϕ×uv1/λ,i=1,2,…,n.After discarding a part of solutions according to the discovery probability, a new solution is generated by(41)xit+1=xit+rxjt−xet,where r is the random number of uniform distribution in (0, 1) interval and xjt and xet are random solutions of the tth generation.


### 5.2.2. Solution of Nonlinear Localization Equations Based on Cuckoo Search Algorithm

The localization equations consisting of equations (36) and (37) are recorded as(42)f1x,y=0,f2x,y=0.When equation (42) is transformed into the unconstrained optimization problem, the objective function of the cuckoo search algorithm is(43)Fw=∑i=12fiw2,where w=x,yT. So solving equation (42) is to find the minimum solution of equation (43).To sum up, the solution process of CS is as follows:(1)
Initialization. The population number of the cuckoo algorithm is m, the maximum number of iterations is K, and the discovered probability is Pa. Moreover, m nest positions are randomly generated and the optimal nest position of the current generation,w0,bestbest∈1,2,…,m, and the optimal fitness value, Fmin, are found. In the simulation experiment of the paper, m = 20, K = 2000, and Pa = 0.25.(2)
Update of Levy Flight Path (Loop Body). The optimal value of the previous generation, wk−1,best, is reserved. K represents the current evolutionary generation. And a new nest location is generated by equation (40). Then we compare the optimal fitness value of the new nest location with that of the previous nest location and keep the solution with the smaller optimal fitness value which is expressed as Ik=wk,1,wk,2,…,wk,m to the next nest location.(3)
Update of the Discovered Probability (Loop Body). For each individual in qk=wk,1,wk,2,⋯,wk,m, r is randomly generated as the probability that the bird’s egg is found by the nest owner, that is, the discovered probability. If r < Pa, the nest value is unchanged; otherwise, the nest is updated according to equation (41). The new nest is compared with qk=wk,1,wk,2,⋯,wk,m, and then the better solution which is expressed as qk=wk,1,wk,2,…,wk,m is kept to the next generation.(4)
Judgment Body. The optimal nest position wk,best in qk and the fitness value Fmin of the solution are calculated. At the same time, it is determined whether the termination condition of the algorithm (the generation is equal to K) is satisfied. If it is satisfied, the global optimal solution is obtained. Otherwise, step 2 is returned, and the iterative operation continues until the termination condition is satisfied.


## 5.2.1. Cuckoo Search Algorithm

The cuckoo search (CS) algorithm is generated by the simulation study of cuckoo’s nesting behavior by Yang Xin-she and DEB Suash. The individuals in birds flocking select nests by Levy flight and abandon by discovered probabilityPa, so the birds population is updated. The paths and locations of cuckoo nests are updated as follows:(38)xit+1=xit+α0xit−xbestt⊕Levyλ,i=1,2,…,n,where xit=xi,1t,…,xi,dtT is the ith solution of the tth generation, α0=0.01, xbestt is the optimal solution of the tth generation, Levyλ∼ϕ×u/v1/λ is the dot product, n is the number of nests, and Levyλ∼ϕ×u/v1/λ is the Levy random search path whose moving step size obeys levy probability distribution as(39)Levyλ∼ϕ×uv1/λ,where u and v follow normal distribution and ϕ=Γ1+λsinλπ/2/Γ1+λ/2×21+λ/2λ, where Γ is the standard Gamma function.Combining equations (38) and (39), in CS algorithm, the new solution xit+1 is generated by(40)xit+1=xit+0.01xit−xbestt⊕ϕ×uv1/λ,i=1,2,…,n.After discarding a part of solutions according to the discovery probability, a new solution is generated by(41)xit+1=xit+rxjt−xet,where r is the random number of uniform distribution in (0, 1) interval and xjt and xet are random solutions of the tth generation.


## 5.2.2. Solution of Nonlinear Localization Equations Based on Cuckoo Search Algorithm

The localization equations consisting of equations (36) and (37) are recorded as(42)f1x,y=0,f2x,y=0.When equation (42) is transformed into the unconstrained optimization problem, the objective function of the cuckoo search algorithm is(43)Fw=∑i=12fiw2,where w=x,yT. So solving equation (42) is to find the minimum solution of equation (43).To sum up, the solution process of CS is as follows:(1)
Initialization. The population number of the cuckoo algorithm is m, the maximum number of iterations is K, and the discovered probability is Pa. Moreover, m nest positions are randomly generated and the optimal nest position of the current generation,w0,bestbest∈1,2,…,m, and the optimal fitness value, Fmin, are found. In the simulation experiment of the paper, m = 20, K = 2000, and Pa = 0.25.(2)
Update of Levy Flight Path (Loop Body). The optimal value of the previous generation, wk−1,best, is reserved. K represents the current evolutionary generation. And a new nest location is generated by equation (40). Then we compare the optimal fitness value of the new nest location with that of the previous nest location and keep the solution with the smaller optimal fitness value which is expressed as Ik=wk,1,wk,2,…,wk,m to the next nest location.(3)
Update of the Discovered Probability (Loop Body). For each individual in qk=wk,1,wk,2,⋯,wk,m, r is randomly generated as the probability that the bird’s egg is found by the nest owner, that is, the discovered probability. If r < Pa, the nest value is unchanged; otherwise, the nest is updated according to equation (41). The new nest is compared with qk=wk,1,wk,2,⋯,wk,m, and then the better solution which is expressed as qk=wk,1,wk,2,…,wk,m is kept to the next generation.(4)
Judgment Body. The optimal nest position wk,best in qk and the fitness value Fmin of the solution are calculated. At the same time, it is determined whether the termination condition of the algorithm (the generation is equal to K) is satisfied. If it is satisfied, the global optimal solution is obtained. Otherwise, step 2 is returned, and the iterative operation continues until the termination condition is satisfied.


## 6. Simulation and Performance Analysis

### 6.1. Simulation Parameter Setting

The American Thunder AN/SPS49 air-to-air search radar is taken as the simulation object in this paper. The main parameters of the radar are shown in Table1.Table 1
Main performance parameters of the radar AN/SPS49.
Parameters
Index
Frequency
f
c
=
850
∼
942
MHz
Transmission power
Peak
P
s
=
360
kW;
Sweep azimuth of radar antenna
Π
=
360
∘
Azimuth rotation rate
Ω
=
36
∘
/
s
3 dB azimuth beamwidth of radar pulse signal
θ
3
=
3.3
∘
Radar pulse repetition frequency
PRF
=
1000
Hz
Rotation angle resolution
Δ
θ
=
Ω
/
PRF
=
0.036
∘
Number of received main beam pulses in one radar sweep cycle
N
p
=
θ
3
/
Δ
θ
=
91.6
Radar pulse length
τ
=
125
μ
sThe horizontal pattern of the antenna main lobe beam is of Gaussian type, and the classical mathematical model as equation (44) is used in the simulation:(44)Gt=exp−4ln2Ωθ32t2.Without loss of generality, the main performance parameters of UAV and receiving loads are set in Table2.Table 2
Main performance parameters of UAV and receiving loads.
Parameters
Index
Receiver sampling frequency
f
s
=
2
MHz
UAV speed
v
=
30
m
/
s
Initial position of UAV
S
i
=
0
,
0
T
Error variance of UAV position measurement
σ
s
2
=
0
m
Main valve pulse receiving interval time
Δ
T
i
,
j
=
60
s
Receiving pulse SNR
SNR
=
20
dBIn this section, the localization accuracy for the radar antenna by the above algorithm is verified. The localization accuracy is calculated byMSES0=∑l=1LS^0l−S02/L, where S^0l is the estimated value of the lth Monte Carlo simulation of the radar antenna localization, L is the number of simulation, and 1000 is selected here.


### 6.2. Sensitivity ofS0 Estimation to Pulse Repetition Frequency and Received SNR

The radar pulse energy received by UAV is inversely proportional to the quadratic ofdi. According to the radar average power, there will also be a higher SNR. Therefore, the radar radiation source is simulated at different positions, and the three target positions are selected as S0,1=60,80Tkm, S0,2=90,90Tkm, and S0,3=160,145Tkm. The type of noise, as discussed in Section 4, is additive random noise with the mean of 0 and the variance of σξ2.Let SNR be monotonically increasing in the interval of [20, 40] dB. When the radar pulse repetition frequency isPRF=1000Hz,PRF=2000Hz, the pulse signal time width is τ=125μs, the UAV speed is v=50m/s, UAV flight time for the single localization ΔTj,i=ΔTj,l=60s, and the localization of three radars are S0,1=60,80Tkm, S0,2=90,90Tkm, and S0,3=160,145Tkm. Three radars are positioned separately. That is, the UAV receives only one radar signal at a time, and there is no signal aliasing. In order to verify the accuracy and effectiveness of the proposed algorithm, Hao’s method in [21] is compared because Hao’s method and the method in this paper both give the solution of the VTDOA localization algorithm. MSE (S0) with different SNRs and PRFs is shown in Figure 4.Figure 4
Estimation accuracy performance as the received SNR increases with different PRFs.
(a)
(b)
(c)As shown in Figure4, with the increase of SNR, the localization accuracy of Hao and the proposed algorithm also become higher. When the noise power becomes smaller, the influence of Δti+1,i becomes smaller, which is more in line with the neglect condition of Δti+1,i. So the localization accuracy is higher. When the PRF is increased, the UAV can obtain more useful information in the same time. So the localization accuracy is higher. As shown in Figure 4, with the increase of PRF, MSE (S0) is reduced and the localization accuracy is higher, which is consistent with the theoretical analysis. By comparing and analyzing the simulation curve of the proposed algorithm and Hao’s algorithm, it can be seen that the localization accuracy is improved under different PRFs and SNRs with the proposed algorithm in this paper. However, for the same target radar, PRF = 1000 Hz and PRF = 2000 Hz, the proposed algorithm has little improvement compared with Hao’s algorithm. This shows that with the increase of PRF, the improvement of accuracy with the proposed algorithm is not obvious compared with that of Hao’s algorithm. This is the result of the estimation of Ψj,i in this algorithm is the same as that in Hao’s. However, with the difference of di, the improvement of Hao’s accuracy will be different which will be analyzed and explained below.When PRF = 2000 Hz, the curves of localization accuracy of three radars versus SNR are shown in Figure5. As shown in Figure 5, when PRF and SNR are the same, the MSE (S0) becomes lower as di increases. Because as di increases, SNB accepted by UAV increases. But the localization accuracy of this algorithm is higher than that of Hao’s algorithm. With different di the improvement of localization accuracy by the proposed algorithm compared with Hao’s algorithm is shown in Figure 6. It can be seen that the improvement becomes better as di decreases. So the proposed algorithm works better in the relatively close localization area, which is also the requirement of the small aircraft to detect the target at close range.Figure 5
Localization accuracy with differentDi.Figure 6
The improvement of localization accuracy with differentDi by the proposed algorithm in this paper.


### 6.3. Sensitivity ofS0 Estimation to li

By changingΔTj,i or v of UAV, the distance li between two UAVs in a single measurement is changed, where li=ΔTj,iv. In the simulation, v=50m/s, which is kept unchanged, li is changed by changing ΔTj,i. The accuracy varies with ΔTj,i at PRF = 1000 Hz, SNR = 30 dB, and τ=125μs as shown in Figure 7. When ΔTj,i increases from 10 s to 150 s, the localization accuracy of the three radars have the same trend. In Hao’s algorithm, when ΔTj,i∈50s,100s, the best localization performance is obtained. When ΔTj,i<50s, with the increase of ΔTj,i, the accuracy is greatly improved. Because the aperture of the VTDOA localization system increases in a reasonable range, it will help to improve accuracy. But when ΔTj,i>100s, the accuracy becomes worse. This is due to that the wavefront of radar signal received by UAV twice will not be the same plane wave because of the excessive li. Therefore, the calculation method of rotation time difference of the radar antenna in Hao’s algorithm brings more errors. In this paper, the NVTDOA localization equations are proposed by using the improved rotation time difference of the radar antenna calculation method. The Newton iteration is used to obtain the convergence solution, which avoids the approximate calculation in Hao’s method. And ΔMSES0 is not less than 100.Figure 7
Localization accuracy withΔT.By the proposed algorithm, with the increase ofΔTj,i, the accuracy increases, and there is no case that the accuracy deteriorates when ΔTj,i is too large. In order to further illustrate the advantages of the proposed algorithm, the improvement of localization accuracy by the algorithm compared with Hao’s algorithm is shown in Figure 8. With ΔTj,i increasing, the accuracy increases, which is consistent with the above analysis. The accuracy with smaller Di is better, so the accuracy of close-range target detection by UAV is better improved in this paper.Figure 8
The improvement of localization accuracy withΔT by the proposed algorithm in this paper.


## 6.1. Simulation Parameter Setting

The American Thunder AN/SPS49 air-to-air search radar is taken as the simulation object in this paper. The main parameters of the radar are shown in Table1.Table 1
Main performance parameters of the radar AN/SPS49.
Parameters
Index
Frequency
f
c
=
850
∼
942
MHz
Transmission power
Peak
P
s
=
360
kW;
Sweep azimuth of radar antenna
Π
=
360
∘
Azimuth rotation rate
Ω
=
36
∘
/
s
3 dB azimuth beamwidth of radar pulse signal
θ
3
=
3.3
∘
Radar pulse repetition frequency
PRF
=
1000
Hz
Rotation angle resolution
Δ
θ
=
Ω
/
PRF
=
0.036
∘
Number of received main beam pulses in one radar sweep cycle
N
p
=
θ
3
/
Δ
θ
=
91.6
Radar pulse length
τ
=
125
μ
sThe horizontal pattern of the antenna main lobe beam is of Gaussian type, and the classical mathematical model as equation (44) is used in the simulation:(44)Gt=exp−4ln2Ωθ32t2.Without loss of generality, the main performance parameters of UAV and receiving loads are set in Table2.Table 2
Main performance parameters of UAV and receiving loads.
Parameters
Index
Receiver sampling frequency
f
s
=
2
MHz
UAV speed
v
=
30
m
/
s
Initial position of UAV
S
i
=
0
,
0
T
Error variance of UAV position measurement
σ
s
2
=
0
m
Main valve pulse receiving interval time
Δ
T
i
,
j
=
60
s
Receiving pulse SNR
SNR
=
20
dBIn this section, the localization accuracy for the radar antenna by the above algorithm is verified. The localization accuracy is calculated byMSES0=∑l=1LS^0l−S02/L, where S^0l is the estimated value of the lth Monte Carlo simulation of the radar antenna localization, L is the number of simulation, and 1000 is selected here.


## 6.2. Sensitivity ofS0 Estimation to Pulse Repetition Frequency and Received SNR

The radar pulse energy received by UAV is inversely proportional to the quadratic ofdi. According to the radar average power, there will also be a higher SNR. Therefore, the radar radiation source is simulated at different positions, and the three target positions are selected as S0,1=60,80Tkm, S0,2=90,90Tkm, and S0,3=160,145Tkm. The type of noise, as discussed in Section 4, is additive random noise with the mean of 0 and the variance of σξ2.Let SNR be monotonically increasing in the interval of [20, 40] dB. When the radar pulse repetition frequency isPRF=1000Hz,PRF=2000Hz, the pulse signal time width is τ=125μs, the UAV speed is v=50m/s, UAV flight time for the single localization ΔTj,i=ΔTj,l=60s, and the localization of three radars are S0,1=60,80Tkm, S0,2=90,90Tkm, and S0,3=160,145Tkm. Three radars are positioned separately. That is, the UAV receives only one radar signal at a time, and there is no signal aliasing. In order to verify the accuracy and effectiveness of the proposed algorithm, Hao’s method in [21] is compared because Hao’s method and the method in this paper both give the solution of the VTDOA localization algorithm. MSE (S0) with different SNRs and PRFs is shown in Figure 4.Figure 4
Estimation accuracy performance as the received SNR increases with different PRFs.
(a)
(b)
(c)As shown in Figure4, with the increase of SNR, the localization accuracy of Hao and the proposed algorithm also become higher. When the noise power becomes smaller, the influence of Δti+1,i becomes smaller, which is more in line with the neglect condition of Δti+1,i. So the localization accuracy is higher. When the PRF is increased, the UAV can obtain more useful information in the same time. So the localization accuracy is higher. As shown in Figure 4, with the increase of PRF, MSE (S0) is reduced and the localization accuracy is higher, which is consistent with the theoretical analysis. By comparing and analyzing the simulation curve of the proposed algorithm and Hao’s algorithm, it can be seen that the localization accuracy is improved under different PRFs and SNRs with the proposed algorithm in this paper. However, for the same target radar, PRF = 1000 Hz and PRF = 2000 Hz, the proposed algorithm has little improvement compared with Hao’s algorithm. This shows that with the increase of PRF, the improvement of accuracy with the proposed algorithm is not obvious compared with that of Hao’s algorithm. This is the result of the estimation of Ψj,i in this algorithm is the same as that in Hao’s. However, with the difference of di, the improvement of Hao’s accuracy will be different which will be analyzed and explained below.When PRF = 2000 Hz, the curves of localization accuracy of three radars versus SNR are shown in Figure5. As shown in Figure 5, when PRF and SNR are the same, the MSE (S0) becomes lower as di increases. Because as di increases, SNB accepted by UAV increases. But the localization accuracy of this algorithm is higher than that of Hao’s algorithm. With different di the improvement of localization accuracy by the proposed algorithm compared with Hao’s algorithm is shown in Figure 6. It can be seen that the improvement becomes better as di decreases. So the proposed algorithm works better in the relatively close localization area, which is also the requirement of the small aircraft to detect the target at close range.Figure 5
Localization accuracy with differentDi.Figure 6
The improvement of localization accuracy with differentDi by the proposed algorithm in this paper.


## 6.3. Sensitivity ofS0 Estimation to li

By changingΔTj,i or v of UAV, the distance li between two UAVs in a single measurement is changed, where li=ΔTj,iv. In the simulation, v=50m/s, which is kept unchanged, li is changed by changing ΔTj,i. The accuracy varies with ΔTj,i at PRF = 1000 Hz, SNR = 30 dB, and τ=125μs as shown in Figure 7. When ΔTj,i increases from 10 s to 150 s, the localization accuracy of the three radars have the same trend. In Hao’s algorithm, when ΔTj,i∈50s,100s, the best localization performance is obtained. When ΔTj,i<50s, with the increase of ΔTj,i, the accuracy is greatly improved. Because the aperture of the VTDOA localization system increases in a reasonable range, it will help to improve accuracy. But when ΔTj,i>100s, the accuracy becomes worse. This is due to that the wavefront of radar signal received by UAV twice will not be the same plane wave because of the excessive li. Therefore, the calculation method of rotation time difference of the radar antenna in Hao’s algorithm brings more errors. In this paper, the NVTDOA localization equations are proposed by using the improved rotation time difference of the radar antenna calculation method. The Newton iteration is used to obtain the convergence solution, which avoids the approximate calculation in Hao’s method. And ΔMSES0 is not less than 100.Figure 7
Localization accuracy withΔT.By the proposed algorithm, with the increase ofΔTj,i, the accuracy increases, and there is no case that the accuracy deteriorates when ΔTj,i is too large. In order to further illustrate the advantages of the proposed algorithm, the improvement of localization accuracy by the algorithm compared with Hao’s algorithm is shown in Figure 8. With ΔTj,i increasing, the accuracy increases, which is consistent with the above analysis. The accuracy with smaller Di is better, so the accuracy of close-range target detection by UAV is better improved in this paper.Figure 8
The improvement of localization accuracy withΔT by the proposed algorithm in this paper.


## 7. Conclusion

When the distance between localization points of UAV is large, the localization error of the VTDOA algorithm increases. First, based on the shortcomings of VTDOA, the NVTDOA localization model of periodic radiation source based on single-station small-load UAV is established considering the integral characteristics of radar antenna rotation time difference. Second, the NVTDOA localization equation is constructed and the convergence solution of the NVTDOA localization equation is obtained by the Cuckoo search algorithm to realize the localization of the radar radiation source. Last, the simulation and comparative analysis of the influence of positioning distance, pulse repetition frequency, SNR, and UAV linear distance on location precision are carried out. The simulation results show that Cuckoo search NVTDOA can improve the positioning accuracy and some conclusions are as follows:(1)
Comparing the localization accuracy of three different position radars, we can see that, in the relatively close localization area, the improvement effect is better. So the NVTDOA algorithm has higher practical value in the close-range target detection by UAV.(2)
With different SNRs, the localization accuracy of NVTDOA is better than that of VTDOA. Therefore, NVTDOA has stronger noise robustness.(3)
The simulation experiments of different motion time of the UAV are carried out. The results show that the smaller the linear distance between two positioning of UAV is, the larger the localization error of the VTDOA algorithm is. But localization accuracy of NVTDOA has little dependence on UAV localization distance, which shows that NVTDOA is more suitable for localization of UAV with high speed and long motion time.



---

*Source: 1410692-2020-03-14.xml*