Backward Projection Positioning for Airborne Moving Targets Using a Distributed Passive Radar System

Abstract

In modern electronic reconnaissance, passive localization plays a crucial role in target positioning and tracking in complex electromagnetic environments due to its inherent concealment capability. However, conventional multistation localization methods often suffer from large estimation errors and limited robustness in dynamic scenarios because of their reliance on multi-snapshot measurements. To address these limitations, this paper extends the backward projection positioning (BPP) method to distributed radar systems for airborne target localization and tracking. By performing backward projection on spatial-domain Doppler history, a cost function involving azimuth angle, elevation angle, and slant range is constructed to generate high-resolution localization images. The distributed radar system forms a physical aperture, enabling high-precision and robust tracking performance. To reduce ambiguities caused by sparse deployment, an optimization strategy based on a genetic algorithm is proposed for radar node placement. The performance of the proposed method is further evaluated using a linear distributed configuration. Simulation results demonstrate that the proposed approach achieves approximately an order-of-magnitude improvement in localization accuracy compared with conventional DOA-based methods.

Graphical abstract

Keywords

passive localization / distributed radar system / backward projection positioning / genetic algorithm

1 Introduction

Traditional active radar is prone to detection and interference, with high costs and limited deployment. With the continuous advancement of electronic systems and the increasingly complex electromagnetic environment on the battlefield, traditional active radar detection faces significant threats and challenges. Passive radar operates by receiving signals without actively transmitting electromagnetic waves. As a result, it has attracted considerable attention due to its superior concealment, longer detection range and higher survivability [1]. It has been widely applied in areas such as radar sensing, navigation, and surveillance [2,3].

On the basis of traditional passive radar, distributed passive radar systems utilize multiple spatially separated receiving stations to detect and localize targets by receiving and processing non-cooperative external electromagnetic signals that are reflected or emitted by the targets [4]. The receiving stations are typically deployed at different geographical locations, enabling target observation from multiple angles and providing richer target information, which facilitates more accurate target identification and localization. Furthermore, the distributed architecture enhances system robustness, allowing continued operation even if some receiving stations are disabled or fail, thereby providing strong resilience against potential damage or system degradation [5].

Traditional non-cooperative radiation source localization methods can be broadly classified into two categories: two-step localization method and direct localization method. The two-step positioning method first involves estimating parameters such as the difference in time of arrival (TOA) [6], frequency of arrival (FOA) [7], direction of arrival (DOA)[8] and their combinations [9,10] of the direct wave signals at each station. Subsequently, the geometric configuration of the distributed passive radar stations is used to construct a corresponding system of positioning equations, and the target location is determined by solving this system. Direct positioning is achieved by fusing the original received signals from all stations at the signal level, followed by estimating the target position based on maximum likelihood principles and spatial grid search techniques [11].

Recently, passive synthetic aperture position (SAP) technology has been proposed by researchers. Wang et al. [12] introduced a synthetic aperture localization model that incorporates the concept of synthetic aperture imaging into passive localization. By utilizing the motion of the platform to form a long virtual aperture, the data are processed in a manner analogous to synthetic aperture radar (SAR), enabling high azimuth localization accuracy. Li et al. [13] proposed a method that utilizes PSA to measure the zenith-passing time of a satellite platform. Yang et al. [14] proposed a multi-subaperture backward projection positioning (BPP) method capable of achieving localization in squint scenes.

Traditional passive localization methods typically require widely distributed receiving stations, and their positioning accuracy is highly dependent on the geometric configuration of these stations. Moreover, these methods generally perform instantaneous localization for moving targets, making them vulnerable to measurement errors, environmental uncertainties, and station failures, which results in limited robustness. This paper extends the BPP method to distributed radar system to achieve high-precision localization of airborne moving targets. The core idea of this work is based on the equivalence between temporal sampling in single-station synthetic aperture systems and spatial sampling in distributed radar systems. This equivalence holds under the following conditions.

First, all distributed radar nodes are assumed to be time- and phase-synchronized to a common reference, ensuring coherent processing of received signals. Second, the target is assumed to be within the common field of view of all radar stations during each observation snapshot, such that spatially distributed measurements correspond to simultaneous observations of the same target state. Third, the propagation environment is assumed to be stationary within each sampling interval, allowing the received phase differences to be interpreted as spatial phase variations analogous to temporal Doppler evolution. This approach enables coherent signal integration, thereby enabling accurate positioning even in complex electromagnetic environments.

The contributions of this work can be summarized as follows:

1) A novel perspective establishing the equivalence between temporal sampling in single-station systems and spatial sampling in distributed systems, thereby providing a theoretical foundation for extending back-projection methods;

2) An adaptation of the back-projection imaging framework to distributed passive radar systems, enabling high-resolution and robust target localization;

3) A node placement optimization method for distributed radar systems based on a genetic algorithm, which effectively suppresses grating lobes and eliminates localization ambiguities.

The remainder of this paper is organized as follows. Section 2 introduces the signal model and presents the BPP method along with the node placement optimization strategy. Section 3 analyzes the performance of the proposed method. Section 4 verifies the effectiveness of the proposed approach through simulations. Section 5 concludes the paper.

2 Localization Model

2.1 Signal Model

A physical long aperture is formed by multiple stations distributed over a wide area. As shown in Fig. 1, the distributed multi-station radar system consists of only N receiving stations deployed on the ground. Due to the large aperture of the distributed radar system, the plane-wave assumption is no longer valid, and the electromagnetic waves must be modeled as spherical waves. The two-dimensional Cartesian coordinates of the receiving stations are denoted as {\mathit{P}}_n=(x_n,y_n,0),n=1,2,\cdots ,M,N. It is assumed that each receiving station shares a common reference time and phase. Furthermore, an airborne moving target \mathit{P}(t)=(x(t),y(t),z(t)) is assumed to be located within the common field of view of all receiving stations, and the stations transmit the captured electromagnetic signals to a fusion processing center. Consider that the airborne moving target transmits a signal as follows

(1)x(t)=s(t)\exp (\mathrm{j}2\text{π}{f}_{c}t)

where s(t) denotes the complex envelope of the transmitted signal x(t) and f_c is the carrier frequency.

The signal received at the n-th receiving station from the target can be expressed as

(2)r_n^{{}^{\prime }}(t)=a_{n}x(t-t_0-{\tau }_n(\mathit{P}))+{\omega }_n^{{}^{\prime }}(t)

where a_n represents the attenuation coefficient of the signal along the propagation path from the target to the n-th receiving station, t_0 denotes the unknown emission time of the signal, {\omega }_n^{{}^{\prime }}(t) denotes additive Gaussian white noise, and {\tau }_n represents the propagation delay from the target to the n-th receiving station, which is expressed as

(3){\tau }_n(\mathit{P})=\frac{\Vert \mathit{P}-{\mathit{P}}_n{\Vert }_2}{c}

where c denotes the propagation speed of electromagnetic waves.

The received signal after down-conversion can be expressed as

(4)\begin{array}{rl}{r}_n^{{}^{″}}(t)=& a_n\exp \{\mathrm{j}\varphi \}s(t-{\tau }_n(\mathit{P}))\times \\ & \exp (-\mathrm{j}2\text{π}{f}_c{\tau }_n(\mathit{P}))+{\omega }_n^{{}^{″}}(t)\end{array}

where \varphi =-2\text{π}{f}_c{t}_0, and {\omega }_n^{{}^{″}}(t) denotes additive Gaussian white noise.

2.2 Backward Projection Positioning

The demodulated signal corresponds to the extracted Doppler component and can be expressed as

(5)r_n(t)=A(n)\exp (-\mathrm{j}2\text{π}{f}_c{\tau }_n(\mathit{P}))+{\omega }_n(t)

where A(n) denotes the complex amplitude, and {\omega }_n(t) represents zero-mean Gaussian white noise with variance {\sigma }^2.

The received signals from all receiving stations can be arranged into a vector as follows

(6)\mathit{R}(t)=\left[\begin{array}{cccc}{r}_1(t)& r_2(t)& \cdots & r_N(t)\end{array}\right]

The BPP method constructs a cost function over the azimuth angle, elevation angle and slant range of the target by leveraging the known positions of the receiving stations. First, the surveillance area is divided into a discrete grid. Then, the Doppler signals corresponding to each grid point are computed. Finally, these computed Doppler signals are coherently combined with the received Doppler signals to generate a localization image of the surveillance region. The detailed procedure is summarized as follows:

1) Grid Partitioning

Grid partitioning is performed on the region of interest, as shown in Fig. 2. Each grid point is parameterized as T=(R_j,{\alpha }_k,{\theta }_l), which corresponds to the Cartesian coordinates (x_{j,k,l},y_{j,k,l},z_{j,k,l}). One radar unit is selected as the reference station, with coordinates (x_0,y_0,0). R_j denotes the slant range from the grid point to the reference radar unit. The azimuth angle {\alpha }_k is defined as the horizontal angle of the target with respect to the reference station, measured from the x-axis in the horizontal plane. The elevation angle {\theta }_l represents the vertical angle of the target relative to the horizontal plane. The corresponding Cartesian coordinates are given by

(7)\begin{array}{l}\{\begin{array}{l}{x}_{j,k,l}=x_0+R_j\cos {\theta }_l\cos {\alpha }_k\\ y_{j,k,l}=y_0+R_j\cos {\theta }_l\sin {\alpha }_k\\ z_{j,k,l}=R_j\sin {\theta }_l\end{array}\end{array}

2) Backward Projection

Given the positions of the receiving stations (x_n,y_n,0), the slant range from each grid point T_{j,k,l}=(x_{j,k,l},y_{j,k,l},z_{j,k,l}) to the n-th receiving station is given by

(8)R_n(j,k,l)=\sqrt{(x_n-x_{j,k,l}{)}^2+(y_n-y_{j,k,l}{)}^2+z_{j,k,l}^2}

The theoretical Doppler signal at the n-th receiving station is expressed as

(9)r_n(j,k,l)=\exp (-\mathrm{j}2\text{π}{f}_c\frac{R_n(j,k,l)}{c})

3) Coherent Accumulation

At each sampling time, the demodulated received signals from all receiving stations are arranged into a snapshot vector as

(10)r_{m{T}_s}(n)=A(n)\exp (-\mathrm{j}2\text{π}{f}_c{\tau }_n(\mathit{P}))+{\omega }_{m{T}_s}(n)

where T_s denotes the sampling period, m is the discrete time index, and {\tau }_n(\mathit{P}) represents the propagation delay from the target position \mathit{P} to the n-th receiving station.

The phase of the rearranged signal is extracted and coherently accumulated with the theoretical Doppler signal r_n(j,k,l) obtained via backward projection. The resulting correlation value is used as the localization outcome for the grid point [R_j,{\alpha }_k,{\theta }_l]. The imaging function is defined as

(11)I_{m{T}_s}(R_j,{\alpha }_k,{\theta }_l)={\left|\sum _{n=1}^N{r}_n^{\ast }(j,k,l)r_{m{T}_s}(n)dn\right|}^2

4) Target Search

The peaks within the generated localization image are detected, and the corresponding X-, Y- and Z-coordinates at each peak are extracted to determine the target positions.

(12)R_{m{T}_s},{\alpha }_{m{T}_s},{\theta }_{m{T}_s}=\arg \underset{R_j,{\alpha }_k,{\theta }_l}{max}{I}_{m{T}_s}(R_j,{\alpha }_k,{\theta }_l)

The motion trajectory of the target is then obtained by fitting the estimated target positions across all sampling times.

2.3 Optimization of Distributed Radar Unit Placement

Due to the structural characteristics of distributed radar systems, such as multiple stations, sparse configurations, and large apertures, image blurring is a common issue in target localization.

To ensure sufficient resolution in the distributed radar system, the system aperture must be properly constrained, which requires fixing the positions of the radar units at both ends of the array. Furthermore, when the spacing between radar units is small, strong mutual coupling effects may severely degrade system performance. Therefore, a minimum-distance constraint must be imposed on the inter-unit spacing to mitigate such effects. The objective of the optimization is to determine the optimal placement of radar units that suppress image blurring while satisfying the above constraints. Specifically, the design goal is to minimize the peak sidelobe level (PSLL). This problem can be formulated as the following constrained optimization problem

(13)

where L_x and L_y denote the aperture lengths along the x-axis and y-axis, respectively, and d_c represents the minimum allowable spacing between radar units.

The algorithm begins with population initialization to generate the initial generation. The population then evolves iteratively through the following steps. First, the fitness of each individual is evaluated. Based on the fitness values, the crossover and mutation probabilities are adaptively adjusted to balance exploration and exploitation. Next, a selection operation is performed, which includes an augmentation–reduction strategy to enhance population diversity while preserving high-quality individuals. Subsequently, generalized crossover and mutation operators are applied sequentially to generate offspring. This iterative evolutionary process continues until the termination criterion is satisfied. The detailed procedure is summarized as follows:

1) Real-Valued Encoding and Initialization of Population

To ensure that the distributed radar system achieves an effective aperture of L_x\times L_y, the positions of the radar units are defined as

(14)\mathit{P}=(x_m,y_n,0),m=1,2,\cdots ,M,n=1,2,\cdots ,N

The spatial coordinates along the x-axis and y-axis are represented as

(15)\mathit{x}=[x_1,x_2,\cdots ,x_M]=[0,d_{x1}+i{d}_{xc},\cdots ,L_x]

(16)\mathit{y}=[y_1,y_2,\cdots ,y_N]=[0,d_{y1}+i{d}_{yc},\cdots ,L_y]

where {\mathit{d}}_x=\left[\begin{array}{c}{d}_{x1},d_{x2},\cdots ,d_{x(M-2)}\end{array}\right] and {\mathit{d}}_y=[d_{y1},d_{y2},\cdots ,d_{y(N-2)}] denote the inter-element spacings along the x- and y-directions, respectively.

The spacing constraints are defined as

(17)

(18)

where d_{xc} and d_{yc} are the minimum allowable spacing between radar units along the x-axis and y-axis, respectively.

In the population initialization stage, M-2 and N-2 random numbers are generated within the interval [0,S{P}_x] and [0,S{P}_y], respectively, and sorted in ascending order. These numbers are used to construct an individual chromosome, which is represented as \mathit{d}={\left[\begin{array}{c}{\mathit{d}}_x,{\mathit{d}}_y\end{array}\right]}^{\mathrm{T}}. The chromosomes of a population consisting of M\times N individuals can be expressed as \mathit{D}=[{\mathit{d}}_1,{\mathit{d}}_2,\cdots ,{\mathit{d}}_{M\times N}].

2) Fitness Function and Preprocessing

The optimization objective is to minimize the PSLL of the localization image in the distributed radar system. The fitness function of an individual with radar unit positions denoted by \mathit{P} can be formulated as

(19)\text{Fitness}(\mathit{P})=\underset{\mathit{P}}{\text{PSLL}}-b

where b=\gamma max(\underset{\mathit{P}}{\text{PSLL}}),0⩽\gamma ⩽1 is the function adjustment parameter.

Since the chromosome \mathit{d} of an individual in the iteration represents a mapping of the radar unit positions \mathit{P}, preprocessing of the population chromosomes is required prior to computing the population fitness. The mapping relationship is expressed as

(20)\mathit{P}=[\mathit{X},\mathit{Y}{]}^{\mathrm{T}}=\mathit{D}+\mathit{C}

(21)\mathit{X}=[d_{x(m,q)}+m{d}_{xc}{]}_{(M-2)\times Q}

(22)\mathit{Y}=[d_{y(m,q)}+m{d}_{yc}{]}_{(M-2)\times Q}

(23)\mathit{C}=[{\mathit{C}}_x,{\mathit{C}}_y{]}^{\mathrm{T}}

where

(24){\mathit{C}}_x=\left[\begin{array}{c}{\mathit{d}}_{xc}\\ 2{\mathit{d}}_{xc}\\ ⋮\\ (M-2){\mathit{d}}_{xc}\end{array}\right],{\mathit{d}}_{xc}={\left[\begin{array}{c}{d}_{xc},d_{xc},\cdots ,d_{xc}\end{array}\right]}_{1\times Q}

(25){\mathit{C}}_y=\left[\begin{array}{c}{\mathit{d}}_{yc}\\ 2{\mathit{d}}_{yc}\\ ⋮\\ (M-2){\mathit{d}}_{yc}\end{array}\right],{\mathit{d}}_{yc}={\left[\begin{array}{c}{d}_{yc},d_{yc},\cdots ,d_{yc}\end{array}\right]}_{1\times Q}

3) Selection Operator Based on Population Enlargement and Subsequent Reduction

Under the discrete real-valued encoding scheme, the population is generated within the theoretical search space [0,SP]. However, during a single execution of the algorithm, the effective search capability is considerably restricted. The crossover operator only recombines existing gene values between parent chromosomes without introducing new genetic information. Although the mutation operator can generate new gene values, its occurrence probability is typically low, and the resulting perturbations are highly stochastic. Consequently, the practical search space is essentially a discrete set formed by the gene values present in the initial population. The optimization outcome is therefore determined by the combinations of these initial gene values rather than the entire theoretical search space. When the population size is small, the algorithm is particularly prone to being trapped in local optima.

To mitigate the tendency of premature convergence and entrapment in local optima, a selection operator based on population expansion and subsequent contraction is introduced. By incorporating controlled perturbations, this operator effectively enlarges the exploration space, improving global search performance. The population is first replicated, after which the individuals in the initial population and the replicated population undergo pairwise combination. The expanded population is constructed as

(26){\mathit{D}}_{\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{p}}=(\mathit{l}-\mathit{r})\circ \mathit{D}+\mathit{r}\circ \overline{\mathit{D}}

where \overline{\mathit{D}} is the replicated population, \mathit{l} is a 2\times (M+N-4)\times Q matrix with all elements equal to 1, \mathit{r} is a 2\times (M+N-4)\times Q matrix whose elements are random numbers distributed in the interval (0,1), and \mathit{A}\circ \mathit{B} represents the Hadamard product of matrices \mathit{A} and \mathit{B}.

By augmenting the population prior to selection, random perturbations within the population are increased, effectively improving population diversity in the search space. Subsequently, a combined selection strategy integrating the roulette wheel method and elitism is employed. Incorporating historical best and elite individuals, N individuals that satisfy the constraints are selected from {\mathit{D}}_{\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{p}} using the roulette wheel method to serve as the parents for the next generation.

4) Generalized Crossover and Mutation

The generalized crossover operator first determines whether crossover should be performed on two genes in the population according to the crossover probability. If the crossover condition is satisfied, two positions on the corresponding chromosomes are randomly selected, and their values are exchanged. Similarly, the generalized mutation operator first determines whether a mutation operation should be applied to a gene based on the mutation probability. If the mutation condition is met, a position on the chromosome is randomly selected, and its value is replaced with a newly generated random number. After performing the generalized crossover and mutation operations, the data on the chromosomes are altered. Therefore, the chromosomes need to be re-sorted to ensure that their elements are arranged in ascending order.

3 Performance Analysis

Under ideal conditions, the measurement information provided by each radar unit of a distributed system at the same time is equivalent to the measurements of the target taken by a single-station system at different times.

To simplify the analysis while preserving the essential geometric characteristics of the problem, the original two-dimensional distributed radar network is reduced to a one-dimensional linear distributed configuration. In the original scenario, the ground stations are deployed at arbitrary positions in the horizontal plane, resulting in a general two-dimensional spatial distribution, which complicates the performance analysis. To reduce analytical complexity, we assume that all ground stations are constrained to lie along a single axis (without loss of generality, the x-axis). Under this assumption, the system is transformed into a linear distributed array, where each station position is expressed as {\mathit{P}}_n=((n-1)d,0,0),n=-N,\cdots ,-1,0,1,\cdots ,N. The element spacing is d.

The received signals at the m-th sampling time can be expressed as

(27)r_{m{T}_s}(n)=a\exp (\mathrm{j}\text{π}{K}_a{n}^2+\mathrm{j}2\text{π}{f}_{d}n+\mathrm{j}\varphi )+{\omega }_{m{T}_s}(n)

where a denotes the normalized signal amplitude, and n=1,2,\cdots ,N represents the index of radar units.

The equivalent Doppler frequency rate K_a is given by

(28)K_a=\frac{d^2{\cos }^2\theta }{\lambda R}

The equivalent Doppler frequency f_d is given by

(29)f_d=-\frac{d\sin \theta }{\lambda }

3.1 Resolution

The grid size is determined by the resolution, and the grid spacing must not be smaller than the resolution value. The squint angle resolution is expressed as

(30){\rho }_{\theta }=\frac{k_{\theta }\lambda }{L\cos \theta }

where k_{\theta }=0.886, and L is the aperture length.

The frequency rate curve can be characterized by the Fresnel function, and the slant range resolution is obtained by

(31){\rho }_r=\frac{K_r\lambda R^2}{L^2{\cos }^2\theta }

3.2 CRLB

The Fisher information matrix associated with the parameter vector \mathit{\beta }={\left[\begin{array}{cc}{K}_a& f_d\end{array}\right]}^{\mathrm{T}} is

(32)\mathit{I}(\mathit{\beta })=\frac{4{\text{π}}^2{a}^2}{{\sigma }^2}\left[\begin{array}{cc}{\displaystyle \frac{N^5}{5}}& 0\\ 0& {\displaystyle \frac{4N^3}{3}}\end{array}\right]

According to the relationship between (K_a, f_d) and (R, \theta ) in Equations (30) and (31), the Fisher information matrix corresponding to the parameter vector \mathbf{\Psi }={\left[\begin{array}{cc}R& \theta \end{array}\right]}^{\mathrm{T}} can be obtained through the Jacobian transformation

(33)\mathit{I}(\mathbf{\Psi })={\mathit{J}}^{\mathrm{T}}\mathit{I}(\mathit{\beta })\mathit{J}

where \mathit{J} denotes the Jacobian matrix defined by

(34)\mathit{J}=\frac{\mathrm{\partial }\beta }{\mathrm{\partial }\mathrm{\Psi }}=\left[\begin{array}{cc}-{\displaystyle \frac{K_a}{R}}& -2K_a\tan \theta \\ 0& -{\displaystyle \frac{d\cos \theta }{\lambda }}\end{array}\right]

The CRLBs for the parameters R and \theta are obtained from the inverse of the Fisher information matrix \mathit{I}(\mathbf{\Psi }), yielding

(35)\{\begin{array}{l}\mathrm{V}\mathrm{a}\mathrm{r}(R)⩾{\displaystyle \frac{{\lambda }^2{R}^4}{4{\text{π}}^2\mathrm{S}\mathrm{N}\mathrm{R}{d}^4{N}^3{\cos }^4\theta }}(\frac{5}{N^2}+{\displaystyle \frac{3d^2{\sin }^2\theta }{R^2}})\\ \mathrm{V}\mathrm{a}\mathrm{r}(\theta )⩾{\displaystyle \frac{3{\lambda }^2}{16{\text{π}}^2\mathrm{S}\mathrm{N}\mathrm{R}{d}^2{N}^3{\cos }^2\theta }}\end{array}

where \mathrm{S}\mathrm{N}\mathrm{R}=a^2/{\sigma }^2 is the signal-to-noise ratio.

It should be noted that the above derivation is valid under the assumption of a sufficiently large number of radar units, for which the asymptotic properties of the estimators are satisfied.

4 Simulation Analysis

To validate the effectiveness of the proposed algorithm and evaluate its performance, simulations are conducted on a linear distributed radar system. Sixteen radar units are uniformly deployed along a 300 m baseline. The first unit is designated as the reference radar, while the remaining units are treated as identical auxiliary radars. A minimum inter-element spacing of 2 m is enforced. In the simulation, the reference radar is located at the origin of the coordinate system. A moving target is initially positioned at (0,10000) m. The target motion is described by a time-varying trajectory \mathit{r}(t)=(x(t),y(t)), where the velocity vector \mathit{v}(t)=(v_x(t),v_y(t)) evolves according to a predefined kinematic model. The target is assumed to transmit a single-frequency signal with a carrier frequency of 3 GHz. The SNR is set to 10 dB. The proposed algorithm is applied to optimize the radar deployment along the baseline and to estimate the target position over time. The target motion is modeled as

(36)x(t)=60t+300\sin (0.05t)\cos (0.1t)

(37)y(t)=10000+80t+200\cos (0.03t)\cos (0.05t)

The localization results are obtained based on the received signals at the initial time instant. First, simulations are conducted for the uniformly distributed radar system. Due to the sparse deployment of the distributed radar units, significant ambiguities are observed in the localization results, resulting in pronounced grating lobes and false target indications in the localization image. The corresponding localization results are shown in Fig. 3.

Then, the radar deployment is optimized to maximize system performance without prior knowledge of the target location. To this end, the entire surveillance region is exhaustively scanned, and the PSLL over the scan is recorded. The maximum PSLL value across the entire region is then adopted as the fitness measure for each individual in the genetic algorithm. This strategy enables the algorithm to identify radar configurations that minimize sidelobe levels over the entire monitored region. The population size is set to N_p=100, the maximum number of iterations is set to N{}_{\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}}=100, the crossover probability is set to P_c=0.9, and the mutation probability is set to P_m=0.05. The convergence behavior of the optimization process is illustrated in Fig. 4, which depicts the evolution of the PSLL across successive generations.

The optimized radar deployment obtained by the proposed algorithm is shown in Fig. 5. The corresponding localization image is presented in Fig. 6. Compared with the uniformly distributed system, the optimized deployment effectively eliminates ambiguities, resulting in an unambiguous localization image.

To further evaluate the performance of the proposed method, a conventional DOA-based localization approach is implemented for comparison. For the conventional DOA-based localization method, a uniform linear array is adopted for each receiving station. The DOA estimation is performed using the MUSIC algorithm. The array spacing is set to half-wavelength to avoid spatial aliasing. For a moving target, each localization instance can be regarded as a single-snapshot measurement at each time instant. Consequently, the DOA-based approach exhibits limited robustness under dynamic scenarios. Identical SNR is used for both methods to ensure a fair comparison. The target is localized at intervals of 0.5 s. Fig. 7 compares the localization results and fitted trajectories of the proposed BPP method and the conventional DOA-based method. It can be observed that the DOA-based results exhibit scattering and larger deviations from the true trajectory, while the proposed method provides a much smoother and more accurate trajectory reconstruction. Fig. 8 further illustrates the positioning error and trajectory fitting error over time. It can be seen that the proposed BPP method consistently achieves significantly lower errors compared with the DOA method.

To evaluate the robustness of the proposed method under different SNR conditions and non-ideal synchronization, Monte Carlo simulations are conducted. The results are shown in Fig. 9(a) and Fig. 9(b). The results indicate that the positioning error decreases as the SNR increases and the phase synchronization error decreases. As the SNR increases, the error gradually converges to a steady-state value.

To further investigate the influence of key system parameters on localization performance, we first evaluate the effect of the number of radar stations while keeping the baseline aperture fixed. The results are shown in Fig. 10(a). It can be observed that positioning error decreases significantly as the number of radar stations increases. Next, we analyze the effect of baseline aperture size while keeping the number of radar stations unchanged, as shown in Fig. 10(b). The results demonstrate that a larger baseline aperture significantly improves positioning accuracy. This is because an increased aperture enhances angular resolution and improves the sensitivity of the received phase with respect to target position.

To demonstrate the capability of three-dimensional target localization, simulations based on a two-dimensional distributed radar configuration are conducted. Due to the high computational complexity associated with node placement optimization in two-dimensional distributed radar systems, the radar nodes in the simulations are arranged in a uniformly distributed configuration with additional random perturbations introduced to reduce periodicity and mitigate ambiguity. Localization results, as shown in Fig. 11, demonstrat that the proposed method can achieve simultaneous estimation of the target’s azimuth, elevation, and slant range.

5 Conclusion

This paper presents an extension of the BPP method to distributed radar systems for airborne moving target localization and tracking. The ground-based distributed radar nodes form a physical long aperture, enabling high-precision and robust target tracking. By exploiting the backward projection of spatial-domain Doppler history, a cost function involving azimuth angle, elevation angle and slant range of the target is constructed to generate a localization image of the target within the surveillance region at each sampling instant. To mitigate ambiguities caused by sparse radar node deployment, a genetic algorithm is introduced to optimize node placement by performing a global search over the surveillance region, ensuring that targets at any location can be accurately resolved without blurring. Simulation results verify the effectiveness of the proposed method, demonstrating approximately an order-of-magnitude improvement over conventional multi-station DOA-based localization method. The proposed approach provides a promising solution for high-precision and robust airborne target surveillance.

It should be noted that the current study does not consider the impact of synchronization errors among distributed radar nodes. In practical scenarios, the large aperture of a distributed radar system imposes significant constraints on radar node placement due to terrain limitations. For example, certain regions may be unsuitable for radar deployment. Future work will focus on improving the robustness of the proposed method against synchronization errors and incorporating more practical constraints into the genetic algorithm, thereby further enhancing its applicability in real-world distributed radar systems.

References

[1] H. Hmam , K. Dogancay . Passive localization of scanning emitters. IEEE Transactions on Aerospace and Electronic Systems, 2010, 46(2): 944–951

[2] A. Napolitano . An interference-tolerant algorithm for wide-band moving source passive localization. IEEE Transactions on Signal Processing, 2020, 68: 3471–3485

[3] Y. Zhu and S. Zhang, “Passive location based on an accurate doppler measurement by single satellite,” in 2017 IEEE Radar Conference (RadarConf), pp. 1424–1427, 2017.

[4] H. Zhu, C. Wang, and L. Wang, “A target detection method of distributed passive radar without direct-path signal,” Electronics, vol. 12, no. 2, p. 433, 2023.

[5] L. Yan , L. Pallotta , G. Giunta , D. Orlando . Robust target localization for multistatic passive radar networks. IEEE Sensors Letters, 2023, 7(7): 1–4

[6] W. Smith , P. Steffes . Time delay techniques for satellite interference location system. IEEE Transactions on Aerospace and Electronic Systems, 1989, 25(2): 224–231

[7] A. Amar , A. J. Weiss . Localization of narrowband radio emitters based on doppler frequency shifts. IEEE Transactions on Signal Processing, 2008, 56(11): 5500–5508

[8] C. Zhang , W. Wang , X. Hong , J. Zhou . A multi-sources doa localization method based on uav cluster systems. IEEE Transactions on Vehicular Technology, 2024, 73(6): 8681–8692

[9] G. Wang , S. Cai , Y. Li , N. Ansari . A bias-reduced nonlinear wls method for tdoa/fdoa-based source localization. IEEE Transactions on Vehicular Technology, 2016, 65(10): 8603–8615

[10] L. Wu , G. Qin , Y. Zou , M. You , B. Chen , D. Chen . Moving target localization in passive distributed mimo radar systems with unknown transmitter positions. Signal Processing, 2026, 239: 110265

[11] O. Bialer , D. Raphaeli , A. J. Weiss . Maximum-likelihood direct position estimation in dense multipath. IEEE Transactions on Vehicular Technology, 2013, 62(5): 2069–2079

[12] Y. Wang , G. Sun , J. Yng , M. Xing , X. Yang , Z. Bao . Passive localization algorithm for radiation source based on long synthetic aperture. Journal of Radars, 2020, 9(1): 185–194

[13] A. Li , H. Huan , R. Tao , L. Zhang . Passive synthetic aperture high-precision radiation source location by single satellite. IEEE Geoscience and Remote Sensing Letters, 2022, 19: 1–5

[14] J. Yang , H. Huan , R. Tao . Multisubaperture backward projection positioning algorithm for radiation sources by single satellite. IEEE Transactions on Aerospace and Electronic Systems, 2025, 61(2): 1267–1282

Rights & permissions

The associate editor coordinating the review of this manuscript was Dr. Junhua Yang