Passive Radar Array Processing with Non-Uniform Linear Arrays for Ground Target ’ s Detection and Localization

The problem of ground target detection with passive radars is considered. The design of an antenna array based on commercial elements is presented, based on a non-uniform linear array optimized according to sidelobe level requirements. Array processing techniques are applied in the cross-ambiguity function domain to exploit integration gain, system resolution and the sparsity of targets in this domain. A modified two-stage detection scheme is described, which is based on a previously-published one by other authors. All of these contributions are validated in a real semiurban scenario, proving the capabilities of detection, the direction of arrival estimation and the tracking of ground targets in the presence of big buildings that generate strong clutter returns. Detection performance is validated through the probability of false alarm and the probability of detection estimation with specified estimation errors.


Introduction
Passive radar technology is being extensively developed to foresee and prevent the multiple threats that modern society faces.In recent years, the improvement in computing capacity and the availability of new technological solutions have increased the interest in Passive Radar Systems (PRSs) as an alternative option to active radars.PRSs present many advantages with respect to active ones: low development, implementation and maintenance costs, small size, low weight, low probability of interception and avoidance of electromagnetic compatibility or environmental impact problems.All of these features are a result of the use of non-cooperative broadcast communications, radar or radionavigation transmitters, as Illuminators of Opportunity (IoOs), rather than a dedicated radar transmitter [1].
A PRS is a multi-static system that allows multiple configurations depending on the number of considered IoOs and receivers, with the objective of detecting targets and estimating their parameters (such as position and velocity).The radar performance is strongly determined by the exploited IoO waveform and the geometry of the radar scenario.The use of IoOs gives rise to complex processing architectures, but thanks to the availability of Field-Programmable Gate Array (FPGA)-based programmable hardware platforms and Graphic Processing Units (GPU), real-time processing is feasible.
Due to the lack of control over the transmitter, two channels are usually used (Figure 1): a reference channel (to acquire the signal transmitted by the IoO) and a surveillance one (to capture the target echoes).To estimate target dynamics information, delay and Doppler-shifted copies of the reference signal are correlated with the surveillance channel in order to generate the Cross-Ambiguity Function (CAF).The received target echoes are extremely faint; therefore, high gain antennas are required in the acquisition stage.If only a single antenna is used in the surveillance channel, the associated narrow beamwidth (due to the inverse relation between antenna gain and beamwidth) leads to a considerable reduction of the coverage area and makes angular discrimination at the detector level unfeasible.
Array processing and digital beamforming techniques have been proposed as a solution for PRS in several works.An FM circular array was considered in [2,3] to obtain one reference and several surveillance beams.In [4,5], a circular antenna array with elements for the VHF-range (150-350 MHz) and the UHF-range (400-700 MHz) is used to exploit alternatively Digital Audio Broadcasting (DAB) or Digital Video Broadcasting-Terrestrial (DVB-T) signals.Two-stage beamforming methods for low complexity Constant False Alarm Rate (CFAR) detection and estimation of the angle of arrival in DVB-T-based PRS have been proposed and tested in [6,7].
In this paper, instead of a Uniform Linear Array (ULA), a spatially-Non-Uniform Linear Array (NULA) for the surveillance channel of a DVB-T passive radar demonstrator, IDEPAR (Improved Detection Techniques for Passive Radars) PRS, is considered.IDEPAR was designed and developed for terrestrial and maritime scenarios, and its capabilities were evaluated in [8][9][10][11][12][13][14].Taking into consideration the IDEPAR design requirement of using Commercial Off-The-Shelf (COTS) components to the greatest possible extent, the array is based on commercial UHF antennas.
The advantage of the NULA is based on improving the radiation beam pattern of a linear aperture of length L with uniformly-spaced array elements [15,16].In the considered architecture, there are different receiving chains, one per single radiating element, to generate the snapshots (vector composed of complex samples associated with the output of each sensor).Snapshot models can be developed in the frequency or time domain [15].In this paper, the spatial filtering is applied in the frequency domain to take advantage of the processing gain of the PRS.Beamforming techniques are applied in a modified two-stage methodology based on [6,7], which is adapted to the NULA and the terrestrial environments: the first step is focused on improving the target detection in the CAF domain, and the second one consists of DoA estimation using a spectrum-based algorithm.
A real bistatic scenario, with the receiver at the rooftop terrace of the Polytechnic School (University of Alcalá), was considered.This scenario was analyzed in [8], and the study is completed with a more accurate estimation of system coverage.In [8], only free space losses were considered.Taking into consideration that the desired targets are on the ground (cars), these estimations were really rough.In this paper, the WinProp (AWE Communications) software is used to obtain more accurate coverage estimations.Results confirm that the designed NULA based on COTS single radiating elements, in combination with the two-stage beamformer, is capable of improving the 2D localization accuracy of terrestrial targets (cars) in semiurban scenarios, providing a low-cost solution for traffic monitoring.
The paper is structured as follows: In Section 2, a general description of the passive radar operating principle and the DVB-T signal as IoO is presented.In Section 3, the problem of antenna design for passive radars is formulated, and an NULA is designed with an optimization process, solved using a genetic algorithm.In Section 4, spatial filtering based on weighting vectors for NULA is described.In Section 5, the modified two-steps spatial filtering technique is detailed.In Section 6, results are presented to validate the study with real data acquired in a semiurban scenario.Finally, in Section 7, the main conclusions are summarized.

General Description and Operating Principle
The basic architecture of passive bistatic radars is depicted in Figure 2: • The reception stage includes the antenna system, the Radio-Frequency (RF) front-ends and the Analog-to-Digital Conversion (ADC) stages for generating the digital reference and surveillance channels.• The processing stage includes the necessary preprocessing to reject interferences and the matched filtering process.The interference can be composed of clutter (interfering radar echoes caused by reflection from objects other than the target), the selected IoO signal present in the surveillance channel (Direct Path Interference (DPI)) and IoO multipath.• Finally, the output of the matched filter is applied to the detection stage in order to decide between target absence or target presence.To perform the matching filtering process, delay and Doppler-shifted copies of the reference signal are correlated with the surveillance channel to generate the Cross-Ambiguity Function (CAF).
For each Coherent Processing Interval (CPI) of duration T int (s), and a sampling frequency f s (Hz), the discrete time expression of the CAF is presented in Equation (1): where: − N CAF is the number of samples in the CPI (N CAF = T int f s ).− m represents the time bin associated with a delay τ m = m/ f s .− p is the frequency bin corresponding to a Doppler-shift f Dop (p) = f s (p/N CAF ).− s re f [n] is the reference signal at the output of the pre-processing stage.− s[n] is the surveillance signal at the output of the pre-processing stage.
In PRS, one of the key issues is the suppression of the DPI, which correlates perfectly with the reference signal.Although it is confined to the zero-Doppler and zero-range bin of the RDS, it produces range and Doppler sidelobes that are several orders of magnitude greater than the echoes that are sought.Actually, the target echo power is usually much lower (e.g., 80-100 dB) than the DPI power.Due to the variability of the signal to be suppressed, adaptive algorithms are used [17,18].The Extensive Cancellation Algorithm (ECA) is widely used [17,19].

Illuminator of Opportunity Selection: DVB-T Signals
Most of the PRS parameters depend on IoO characteristics: availability, radiated power, coverage area, transmitted waveform.In the passive radar literature, digital commercial waveforms have been extensively analyzed to determine their impact on system resolution, detection ambiguities and expected coverage: DAB [20], Global System for Mobile (GSM) [21], Universal Mobile Telecommunication System (UMTS) [22] and DVB-T, the European Union standard, [23,24].
The Ambiguity Function (AF) is a powerful tool to analyze the detection capabilities of radar signals [25].The AF allows the estimation of the bistatic range and Doppler shift of the target and provides the necessary signal processing gain to allow the target echo detection, performing as a matched filter for the radar system.DVB-T is the digital terrestrial television broadcast transmission standard most used in the European Union [26] and has been used as IoO in this study.This system transmits compressed digital audio, digital video and other data in an MPEG transport stream, using Coded Orthogonal Frequency-Division Multiplexing (COFDM) modulation.COFDM works by splitting the digital data stream into a large number of digital streams, which digitally modulates a set of closely-spaced adjacent sub-carrier frequencies.In the case of DVB-T, there are two options for the number of carriers, known as 2 K-mode or 8 K-mode.
DVB-T signals are composed of a random component, generated by the interleaving and the OFDM modulation of the MPEG-2 signal, and a deterministic component generated by the presence of signaling elements (Guard Interval (GI), Pilot Carriers (PC) and Transport Parameter Signaling (TPS)).Deterministic components of the DVB-T signal give rise to ambiguity peaks in the AF.The range of ambiguities is summarized in Table 1.In the scenario under analysis, located in Madrid (Spain), the 8 K mode was adopted.In Figure 3, the zero Doppler cut of the AF of a single channel DVB-T signal is presented.The expected ambiguity peaks are shown together with the main peak located at the range-Doppler plane origin.

Problem Formulation
The design of an array antenna to be part of the IDEPAR demonstrator [8], which has been used to test the algorithms, is one of the objectives of this paper.Figure 4 presents a functional block diagram of this multichannel DVB-T-based PRS.To ensure an easy maintenance and low-cost, COTS components were used as much as possible, and following this design requirement, commercial antennas were selected.The use of an antenna array and array digital signal processing could provide further benefits: • Reducing the direct signal interference, by creating suitable nulls in the beam pattern in the direction of transmitters.
• Enabling the application of beamforming techniques for estimating the target azimuth location.
• Increasing the angular coverage using an individual element of the array with a wide beamwidth.
The Televés 4G NOVA antenna (Figure 5), which is a log-periodic antenna designed in microstrip technology for operating in the 698 MHz-960 MHz and the 1700 MHz-2700 MHz bands, was selected according to the following requirements: (a) broad beamwidth to provide the required angular coverage and a wide range of steering angles; (b) high discrimination between the main and the back lobes; (c) small physical dimension.The 3D antenna model was generated and simulated in the Ansoft HFSS (High Frequency Structure Simulator). Figure 6 shows the radiation patterns in elevation and azimuth.The 3-dB azimuth beamwidth is approximately 60 • , and the directivity is 7.98 dBi for 770 MHz (Figure 7).In a first approach, a study of Uniform-Linear Arrays (ULA) along the y-axis was carried out.Arrays of five elements were considered because the current version of the passive radar demonstrator includes six synchronized receiving chains: one for the reference channel and five for the surveillance one.In Spain, DVB-T is transmitted in horizontal polarization, which requires all antennas to be parallel to the ground.The inter-element spacing is limited by the width of the antenna, which could give rise to the appearance of grating lobes.Equation ( 2) is the stringent condition for no grating lobes [27], being φ sa the steering angle and d the inter-element spacing.
For d = 315 mm, the antenna width, beam patterns generated for steering angles higher than 13.7 • and lower than −13.7 • will present grating lobes at the operating frequency of 770 MHz.Although the radiation pattern of the antenna reduces the effect of grating lobes, they can be seen in Figure 8 for steering angles greater than or equal to 15 • .A further study revealed that the reflector shown in Figure 5 only worked at the higher frequency band (1700 MHz-2700 MHz), so we decided to remove it together with the radome, reducing the minimum inter-element spacing to 210 mm.In this case, the steering angles' interval for no grating lobes is |φ sa | 58.79 • .
As the single radiating element beamwidth is 60 • , the realistic steering angle interval is |φ sa | 30 • .Therefore, the inter-element spacing can be increased to improve array directivity while fulfilling the no grating lobes condition expressed in Equation (2).If this condition is relaxed allowing the appearance of part of the grating lobe for a steering angle φ sa = 30 • , with a maximum level equal to the Sidelobe Level (SLL), the inter-element spacing can be increased to 270 mm (Figure 12b).In this work, when the SLL and Grating Lobe Level (GLL) are expressed in (dB), they are calculated as the difference in dB relative to the main lobe level.Therefore, the higher the value of SLL or GLL in dB, the lower the importance of these lobes in the radiation pattern.
Simulation results for different steering angles (φ sa = [−30 • , −15 • , 0 • , 15 • , 30 • ]) for ULAs with inter-element distances of 210 mm and 270 mm are summarized in Table 2.The SLL values obtained for φ sa = 30 • with the ULAs are clearly bad.In order to obtain a compromise solution between directivity and SLL, maintaining the beam steering capability in an azimuth range of ±30 • , a solution based on Non-Uniform Linear Arrays (NULAs) was studied.

Design of a Non-Uniform Array for Optimizing the SLL Value
In this subsection, the design of the NULA for the considered five-element array antenna is described.The objective was to minimize the importance of the sidelobe level, so that the absolute value of the difference ( 1 norm) between the normalized radiation pattern and a specific value of sidelobe level was calculated for φ outside the main beam width (BW 3dB ), which defines the system azimuth resolution, and only for the values of φ where the normalized radiation pattern is higher than the specific value of sidelobe level.When the sidelobe level is minimized, BW 3dB increases, and the resolution decreases.In order to control BW 3dB , an additional term was added to the cost function, to get a trade-off between sidelobe level and main beam width.The preliminary cost function to be minimized is presented in Equation (3): • d and k are the parameters to be adjusted: , for an array with N elements), and has been adjusted using a genetic algorithm; k is the parameter that allows us to control the main beam width, BW 3dB , and has been adjusted experimentally.
is the magnitude of the normalized radiation pattern of the array along the y-axis expressed in natural units, particularized in θ = 90 • , the azimuth plane, for the considered distances vector d.Note that for the evaluation of C(d, k, sll max ), only azimuth values outside the main beam and with an associated radiation pattern magnitude higher than sll max are considered in the summation.
• sll max is the specified maximum value of the sidelobe level in natural units.SLL max (dB) is defined as the difference between the main beam level in decibels and the highest sidelobe level in decibels: is the main beam width for the considered distances vector d.
As the array must be pointed towards a set of azimuth directions ranging from −φ max to φ max , a vector of N φ steering angles was defined, and a radiation pattern was generated for each pointing direction using conventional beamforming weights, Ēφ sa,i (φ, d; φ sa,i ).The final cost function was defined as the average of the preliminary one, throughout the azimuth sector of pointing angles defined by Φ sa , Equation (4): Once the cost function was defined, the optimization problem to be solved for an array of N elements was formulated with Equation ( 5), where d min and d max are the limits of the allowed distances between elements.
The optimization problem was solved using a genetic algorithm [28], using MATLAB R , with the default options: • Initial population size: 20 individuals.
• Number of best individuals that survive to next generation without any change: two.
• The fraction of genes swapped between individuals: 0.8.
• The number of generations between the migration of the fittest individuals to other sub-populations: 20.• Fraction of those individuals scoring the best that will migrate: 0.2.
• Number of generations to be simulated: 100.
• Stopping criteria based on the maximum number of generations; if after this number of generations there is no improvement, the simulation will end.
The maximum SLL value was defined to be SLL max = 14 dB.A set of pointing directions Φ sa = {−30 • , −15 • , 0 • , 15 • , 30 • }, was defined.Distances should be greater than d min = 210 mm (the minimum value defined by the antenna width after removing the radome and the reflector) and lower than d max = 270 mm, in order to avoid higher grating lobes in the desired azimuth coverage.
The optimization procedure was applied to estimate the optimal inter-element distances of a five-element array antenna, for different k values.Good results were obtained for 0.   ] and are presented in Table 2.The worst results of SLL (dB) and GLL (dB) are associated with steering angles of 30 • and −30 • .For these angles, the obtained SLL (dB) is lower than SLL max (dB) due to the condition imposed in the optimization cost function related to maintaining a narrow beamwidth, controlled with k = 1.Results confirm that the designed NULA is able to provide the desired azimuth coverage with GLL (dB) and SLL (dB) higher than ULA with d = 270 mm, and maintaining the directivity higher than ULA with d = 210 mm (Figure 11).We have repeated the procedure with an array of 11 elements (10 inter-element distances).The results are presented in Table 3 and in Figure 12.In this case, the resultant distances (d = [349, 259, 223, 210, 210, 210, 214, 235, 262, 231] mm) correspond to an asymmetric NULA that provides a directivity 0.6 dB lower than that associated with an ULA of 11 elements with d = 270 mm, but with an SLL (dB) that is approximately 4 dB higher, which demonstrates the usefulness of the proposed method.

Array Processing
Beamforming can be applied with weight vectors if the narrowband constraint expressed in Equation ( 6) is verified, B being the signal bandwidth, da max the maximum physical distance between two elements in the array and TBWP the time-bandwidth product [15].
Snapshots for array signal processing can be generated in the time domain (outputs of the analog-to-digital converters associated with each sensor) or in the frequency domain (outputs of the CAF processing stage calculated for each sensor).In PRSs, the CAF implements a time-frequency mapping with a time (bistatic delay) resolution that depends on the signal bandwidth and a frequency (Doppler shift) resolution that depends on the integration time.If beamforming is applied in this transformed domain, we could take advantage of the integration gain and the range-Doppler mapping of the targets and interferences present in the input signal [6].Interferences are expected to concentrate around the zero Doppler line of the CAF, while targets will concentrate in the cells associated with their position and speed (range, Doppler), respectively, giving rise to sparse matrices from the point of view of targets.Therefore, if snapshots are defined in the CAF domain, a specific beam space could be generated for each range-Doppler (time-frequency) cell.The CAF will separate targets according to range and Doppler, and in each range-Doppler cell, each beam space will be only responsible for the azimuth discrimination.
Figure 13 shows the processing chain previous to the array signal processing stage.If the snapshots were defined in the time domain, the resulting vector would be: In this paper, we rely on the definition of the snapshots in the frequency domain, and an independent CAF must be calculated from each s i [n], i = 1, ..., N. The snapshots are obtained for each (m, p) in the CAF domain: The processing scheme requires the division of the acquisition time into Coherent Processing Intervals (CPIs), of duration T int .A set of N CAFs is generated for each CPI (Figure 13).The integration time, T int , must be significantly greater than the propagation time across the array T int > da max c [15].The second requirement on CPI is determined by the bandwidth of the input signals and the shape of their temporal spectra.When the product B • T int is large enough and the input signals' temporal spectra are almost flat, the frequency-domain snapshot model is usually applied [15].
The steering vector of a signal impinging from a direction φ sa,i is defined in Equation ( 7), where f s is the sampling frequency, f p is the carrier frequency, N CAF is the number of samples in a CPI and c is the speed of light.

Data Independent Beamforming
Data independent beamforming is usually based on the maximization of array directivity or the minimization of the squared error between the beamformer output and the desired radiation pattern at defined directions for sidelobes control [15,29,30].In this subsection, both approaches are summarized to provide the theoretical basis of the two-stage array signal processing scheme for improving detection and tracking in PRSs, which is presented in Section 5.
• Maximization of array directivity: Taking into consideration the snapshots defined in the CAF domain, array directivity is calculated as where matrix A(p) is calculated in Equation (9).The element at the m-th row, and n-th column can be expressed with Equation (10).
The vector w(p, φ sa,i ), which maximizes D(p, φ sa,i ), is expressed in Equation (11).In the special case of ULA, A(p) tends to the identity matrix, and the maximum directivity weight vector is close to the uniform weight vector steered to the desired direction.For the special case of a ULA with inter-element spacing equal to λ/2, A(p) is equal to the identity matrix.
• Sidelobes reduction: If lower sidelobes are desired, sacrificing some directivity, a methodology based on [31] can be developed.The φ space is divided into R δ sectors, The objective is to find a single weight vector, w SLL , to approximate the radiation pattern obtained as the concatenation of the set of defined target radiation patterns | Ēt,i (φ)|, φ ∈ δ i , i = 1, ..., R δ .This is equivalent to maximize the directivity controlling the SLL level in the whole φ interval, (0, 2π]. If the radiation pattern obtained using w SLL is | ĒSLL (φ)| = w H SLL a(φ), φ ∈ (0, 2π], the squared error to be minimized is defined in Equation ( 12), being If snapshots are defined in the CAF domain, the methodology described so far must be applied for every Doppler shift.Assuming that the target radiation pattern is the same for all Doppler shifts, a weight vector must be determined for each p value.On the other hand, multiple beams pointing towards a set of azimuth directions Φ = {φ sa,1 , φ sa,2 , ..., φ sa,Nφ sa } are required, so the final error function to be minimized is a function of p and φ sa,i as expressed in Equation (13), with: The minimization problem to be solved is formulated in Equation ( 14) where ε 2 max is the imposed SLL constraint.Unfortunately, there is no closed-form solution, and an iterative procedure is required.In this paper, the solution proposed in [31] was implemented.min

Directivity Maximization Applied to ULA and NULA
In Section 3.2, a comparative study of the ULAs and NULAs composed of five and eleven commercial antennas was carried out.Conventional uniform weighting was applied to analyze SLL, GLL and directivity in order to prove the benefits of designing NULAs.In this section, a parallel study is presented applying the maximum directivity weights defined in Equation (11).ULAs' maximum directivity weights are close to the conventional uniform weights, but in the case of NULAs, they are different.Tables 4 and 5 summarize the radiation pattern parameters for the two considered ULAs (with inter-element distances equal to 210 mm and 270 mm) and the NULA designed in Section 3.2 with beamformer weights calculated for maximizing the directivity, for N = 5 and N = 11 elements, respectively.The values obtained for the ULAs are very similar to those presented in Tables 2 and 3, as expected.Even though for N = 5, SLL and directivity of the ULA with d = 270 mm and the NULA are very similar, for N = 11, the SLL provided by the NULA is significantly better, maintaining a similar directivity.In this case, weight vectors are calculated for minimizing Equation (12), and a comparative study of the ULAs and NULAs composed of five and eleven commercial antennas is presented.The objective is to evaluate the impact of this beamforming technique on radiation pattern parameters in both array architectures.The iteration method proposed in [31] was applied with the constraints of SLL and GLL to be lower than 15 dB.Tables 6 and 7 summarize the results for N = 5 and N = 11, respectively.The SLL improvement with respect to uniform conventional beamforming is clear, with a very small reduction in directivity.For N = 5, the SLL (dB) of the NULA is slightly higher than the SLL (dB) of the ULA with d = 270 mm, and the directivity is slightly lower, but GLL (dB) is lower than the SLL (dB) for the ULA with d = 270 mm for steering angles equal to ±30 • and does not fulfill the defined GLL constraint.The gain and SLL (dB) are similar for NULA and the ULA with d = 210 mm, but the BW in the azimuth plane is 2 • narrower in the case of NULA.For N = 11, the improvement of NULA compared to ULAs is significantly higher.The SLL (dB) of the NULA is higher, especially for steering angles close to the broadside; the directivity is lower than the directivity of the ULA with d = 270 mm, but this ULA does not fulfill GLL constraint for steering angles equal to ±30 • .

Two-Stage Spatial Filtering for a PRS
In order to provide 2D target localization (range and azimuth), a two-stage spatial beamforming approach based on CFAR detection in multiple simultaneous beams and on the estimation of the angle of arrival is proposed.Figure 14 depicts the scheme of this processing method.It is based on modifications of those proposed in [6,7] for ULAs and implemented for the NULA designed in Section 3.2: -In [6], Spatial Adaptive Processing (SAP) was analyzed for PRSs, for interference spatial filtering and the estimation of the target's Direction of Arrival (DoA).Time-domain SAP and range-Doppler domain SAP were studied, and an alternative PRS signal processing chain, using both spatial filtering and high resolution DoA estimators working in the range-Doppler domain, was proposed.-In [7], real-time detection and tracking were demonstrated using a two-stage beamforming approach.A smaller set of beams was used in the first stage for CFAR detection.For each [m,p] pair, the beam that provided the maximum power was determined, and the CFAR detector was applied on this beam's outputs.First stage beams were designed under the condition of the maximally flat response between beams' directions.A more crowded set of beams was used in the second stage for DoA estimation using a spatial smoothing-based Bucci algorithm proposed in [32].
The modifications proposed in the two-stage spatial filtering scheme are related to the generation of orthogonal beams using weight vectors calculated under SLL requirements (w SLL ) for CFAR detection, the CFAR detector threshold estimation in the 3D space defined by range, Doppler and azimuth and the application of DoA techniques to estimate the azimuth of the detected targets using weight vectors calculated to maximize the directivity (w MD ): • The first stage implements a beamforming processing using orthogonal beams calculated under SLL requirements.If orthogonal beams are used, a signal arriving along the maximum radiation axis of a beam will have no output in any other beam; and a signal that is not along the maximum radiation axis will appear in the sidelobes of the other beams [15].The set of steering directions Φ SLL = {φ SLL,1 , φ SLL,2 , ..., φ SLL,Nφ SLL } is defined following an iterative process that starts with the broadside beam and continues defining the maximum radiation axis of the side main beams along the nulls of the broadside beam.This process is iterated until the azimuth sector defined by the single radiating element of the NULA is covered.For each φ SLL,i pointing direction and each Doppler shift, p, the optimization problem formulated in Equation ( 14) must be solved to obtain the weight vector w SLL (p, φ SLL,i ), generate the associated beam and determine the main radiation axis of the side main beams.The result is a three-dimensional matrix called S CAF , with indexes m, p and φ SLL,i , for which elements s CAF [m, p, φ SLL,i ] are obtained with Equation ( 15): The CFAR detection is implemented as follows: 1.For each [m, p] pair, the Cell Under Test (CUT) is the cell s CAF [m, p, φ SLL,i ] for the φ SLL,i value where the maximum power of the echo is received.2. A 3D reference window with dimensions R R xR D xN φ SLL (range xDoppler x steering angle), is generated using the neighbor cells around the CUT and excluding the guard cells, which are defined along the range and Doppler for all of the steering angles, to reduce the impact of extended targets in the threshold estimation.This 3D reference window also excludes The use of 3D reference windows instead of 2D or 1D improves interference statistics estimation, reducing the CFAR losses [33].The outputs of this stage are the detected targets and the estimation of their ranges and Doppler values.
A new set of steering angles is defined, Φ DoA = {φ DoA,1 , φ DoA,2 , ..., φ DoA,N φ DoA }, with N φ DoA > N φ sa to increase azimuth estimation accuracy.For each DoA steering direction and Doppler shift, p, the weight vector is calculated to maximize the directivity using Equation (11).DoA estimation is based on the beamformer output spectrum expressed in Equatino (16), for (m, p) pairs where a target was declared.
R s s s s (m, p) is the instantaneous spatial covariance matrix estimation for the snapshot s s (m, p).For each (m, p) pair, the value of φ where the beamformer spectrum, S m,p (φ), is maximum determines the estimated azimuth for the target detected at (m, p) delay and Doppler bins.
As an alternative, a parametric estimator that exploits the information of R s s s s (m, p) eigenvalues and eigenvectors, such as MUSIC (Multiple Signal Classification), could have been used.However, MUSIC resolution decreases when it is applied to NULAs [15].

Scenario Description
The capabilities of the designed NULA and the two-step beamforming processing are evaluated in a real radar scenario, in the nearby area of the Polytechnic School of University of Alcalá (Spain).The PRS is located at the rooftop terrace of the Polytechnic School.This scenario is characterized as a semiurban environment: low-height buildings belonging to the University Campus surrounded by countryside and several roads (Alcalá-Meco road, R2 highway and inner-campus secondary roads).A controlled car with a GPS device was running along Alcalá-Meco road within the area of interest.
In the considered radar scenario, there is one DVB-T transmitter (Torrespaña), with an Equivalent Radiated Power (ERP) equal to 20 kW, which is used as IoO taking into consideration its high radiated power and its omnidirectional radiation pattern.The estimated received power levels provided by this IoO in the surveillance area are depicted in Figure 15.Coordinates of the main elements of the geometry are presented in Table 8.In order to obtain the system coverage, the Bistatic Radar Cross-Section (BRCS) of the targets of interest was estimated in [8].A Mazda 6 SPORT vehicle was considered as a representative target for traffic monitoring.A study of the scenario geometry revealed that the bistatic angles varied from 105 • -135 • .In Table 9, BRCS representative values are summarized.For coverage studies, σ bistatic = 10.6118dBsm was used.In Figure 16, the single element antenna radiation pattern in the horizontal plane is superimposed on the scenario image to show the instrumented coverage sector.The antenna orientation was selected to ensure that the maximum of the radiation pattern was pointed towards the area of interest in order to enclose the roads under study.

Array Signal Processing Performance Evaluation
The NULA with five elements designed with the method described in Section 3.2 fulfills the narrowband condition (TBWP = 0.0261 < 1), and the product B • CPI = 2 × 10 6 is large enough to implement the frequency-domain snapshot model.The two-stage processing proposed in this paper was applied to the designed NULA.The coverage for only one Televés 4G NOVA antenna in the surveillance channel is depicted in Figure 18a, for P D = 0.8 and P FA = 10 −6 , using system parameters summarized in [8].In [8], only basic propagation losses were considered; in this new study, WinProp (AWE Communications) software was used, including geographical data and taking into consideration that the desired targets were located on the ground.Results presented in Figure 18b confirm the range coverage improvement obtained using an array antenna in the surveillance channel.The P FA requirement was set to 10 −5 to determine the adaptive CFAR threshold.The superposition of the detection matrices obtained from the processing of 80 consecutive CPIs is presented in Figure 19.Results provided by only a single radiating element of the array are shown in Figure 19a.Two main groups of detections are identified: a set of trajectories located in the Alcalá-Meco road (marked with the green ellipse); a second set associated with R2 highway (marked with the red ellipse).From visual inspection, the number of detected targets and the detection densities of the guessed trajectories are clearly higher for the NULA and the detection stage array processing.For performing a more rigorous detection performance evaluation, three trajectories were selected without loss of generality (Figure 20): • Trajectory 1 (T1): a cooperative moving target provided with a GPS receiver, running along Alcalá-Meco road, observed from 25th CPI-49th CPI.• Trajectory 2 (T2): a moving target in Alcalá-Meco road, observed from 28th CPI-79th CPI.
• Trajectory 3 (T3): a moving target in R2 highway, observed during the whole acquisition interval (from 0th CPI-79th CPI).To estimate P FA and P D , ground-truths at the output of the detector are required, but due to the complex nature of the electromagnetic propagation processes, targets dynamics and radar system, the real ground-truth is not always available.Using the methodology described in [8], a ground-truth was generated for each CFAR detector using GPS data of cooperative vehicles and visual information about non-cooperative targets present on R-2 highway and Meco road during the acquisitions.The P D was calculated at the plot level due to the errors associated with the pixel detection grouping techniques; P FA was estimated throughout CAF areas where no targets were expected.Monte Carlo techniques were applied, guaranteeing an estimation error lower than 10%.
Table 10 details the estimated P D and P FA for the trajectory associated with the cooperative target, T1.Table 11 presents the number of detections at the plot level for the three targets, T1, T2 and T3.In both cases, results are provided for beamforming weights calculated for controlling SLL, w SLL , and for beamforming weights calculated for maximizing the directivity, w MD .In all cases, the P D and the ratio of detections per target are higher for the weights calculated for controlling the SLL.
Results presented in Tables 10 and 11 require a further explanation because a higher P D was expected for the controlled car, T1, and the number of detections for T1 and T2 was expected to be higher than for T3 (actually, T3 is significantly further than T1 and T2).The main reason responsible for these results is shown in Figure 20: a big building made of aluminum (IMMPA) is affecting the propagation of the electromagnetic waves.Although cars can be seen from the passive radar emplacement, the diffraction losses generated by the IMMPA are really high.This effect can be observed in Table 12, where Signal-to-Interference Ratios (SIR) were estimated using the methodology described in [8].The estimated SIR for targets T1 and T2 is significantly lower than that estimated for target T3, although T3 is further.The (m, p) coordinates in the CAF domain of each declared target are supplied to the second array signal processing stage to estimate their azimuths.The DoA estimation is carried out in a new beam-space with better accuracy (steering angle increment of 0.01 • ) and better resolution, using beamforming weights calculated for maximizing the directivity.Figure 21 depicts the radiation patterns with a maximum radiating axis towards the broadside, for the NULA with weights calculated to maximize the directivity, and weights calculated to control the SLL, showing a narrower main beamwidth for the first one.For the (m, p) pair associated with each declared target, a new set of radiation patterns is generated with maximum radiating angles varying from −30 • -−30 • , with a step of 0.01 • (Φ DoA = {−30 • , −29.99 • ...., 29.99 • , 30 • }); the target DoA is estimated as the maximum of the beamformer output spectrum expressed in Equation (16).
As an example, the detection matrix of the 30th CPI is presented in Figure 22a.The detections associated with each considered target (T1, T2 and T3) are marked in yellow.The beamformer output spectrum built for each target detection is shown in Figure 22b.The estimated target azimuth will be the maximum of the associated beamformer output spectrum.Figure 23 presents the estimated angles of all of the detections for the three considered trajectories in the acquisition interval.At this point, 3D targets parameter space can be transformed to 2D map coordinates.When the trajectories are depicted in the radar scenario (Figure 24), the target movement dynamics fits well with expected target echoes generated by vehicles running along the Alcalá-Meco road and R2 highway, respectively.GPS data are also represented validating the results obtained with the NULA and the two-stage beamformer.

Conclusions
This paper tackles the problem of radar detection in a 3D target parameter space using a DVB-T-based PRS, with an NULA in the surveillance channel and spatial filtering in the frequency domain.The proposed solution allows the estimation of target azimuth, Doppler and range, to improve the target localization accuracy and to make traffic monitoring easier.Independent acquisition chains were considered for digitizing the signal acquired by each single element of the array.Taking into consideration the requirement of using COTS components in the demonstrator to the greatest possible extent, the commercial Televés 4G NOVA antenna was selected as the array element for the surveillance channel.This antenna provides a good range of steering, reduced dimensions, high front-to-back ratio and the capability to fulfill the high gain requirement for passive radar detection.
The use of non-uniform linear arrays was proposed, and a cost function was designed to determine the inter-element distances, to obtain a compromise solution of sidelobe levels and main beam width.The formulated optimization problem was solved using a genetic algorithm.Arrays with five and eleven elements were designed.N = 5 was selected because this is the number of available acquisition chains for the surveillance channel of the passive radar demonstrator.N = 11 was selected to prove the potential improvement achievable when more acquisition chains are available, as expected in operative systems.
Beamforming has been applied in the frequency domain, using a modified version of the previous two-stage algorithms published by other authors: • The first step generates a full-dimensional beam-space based on orthogonal beams in the azimuth coverage area where the CFAR detector is applied.Beamforming weights are calculated for guaranteeing an SLL (dB) lower than a specified value (15 dB in the presented results).In the CFAR detector, the windowing technique defines CUTs related to maxima in the beam-space and 3D reference windows that exclude a set of guard cells to reduce the impact of targets contributions around the CUT, in detection threshold estimation.The use of 3D windows reduces CFAR losses.• The angles of arrival of the detected targets are estimated using the output spectrum of a new denser beam-space generated using weight vectors calculated to maximize the directivity.
Performance improvement associated with the use of NULA instead of ULA and the combination of weighting techniques to control SLL and maximizing directivity, instead of uniform conventional beamforming, was proven by simulation.The designed NULA constitutes a good base architecture to implement beamformers under the specified design criteria, outperforming the ULA for the whole system coverage sector, especially when eleven elements are considered.
Finally, a semiurban radar scenario, with the PRS located at the rooftop terrace of the Polytechnic School of the University of Alcalá, was analyzed to validate the capabilities of the considered solution for monitoring terrestrial targets.The environment is characterized by low-height buildings and several roads.PRS coverage was estimated using WinProp (AWE Communications) software and including geographical data, considering that the desired targets are on the ground.A controlled car with a GPS device was used during the measurement campaign.Results obtained with real data validate the conclusions extracted from the simulations of the proposed array architecture and the two-stage array signal processing.

Figure 1 .
Figure 1.Basic operating scheme of a Passive Radar System (PRS).

Figure 3 .
Figure 3. Zero-Doppler cut of the Ambiguity Function (AF) of a simulated Digital Video Broadcasting-Terrestrial (DVB-T) signal.

Figure 4 .
Figure 4. Functional block diagram of Improved Detection Techniques for Passive Radars (IDEPAR).
(a) Elevation radiation pattern (b) Azimuth radiation pattern

Figure 8 .
Figure 8. Azimuth plane radiation pattern of the ULA with five Televes 4G NOVA antennas, for 770 MHz, inter-element spacing equal to 315 mm and three steering directions: 0 • (broadside), 15 • and 30 • .The effect of grating lobes is marked with an ellipse.

Figure 9 .
Figure 9. Designed NULA based on Televés 4G NOVA antennas without the radome detailing the inter-element spacing.

Figure 11 .
Figure 11.Azimuth plane radiation patterns of the ULAs (d = 210 mm and 270 mm) and NULA, with five Televes 4G NOVA antennas without the radome or reflector, for 770 MHz.

Figure 12 .
Figure 12.Azimuth plane radiation patterns of the ULAs (d = 210 mm and 270 mm) and NULA, with 11 Televes 4G NOVA antennas without the radome or reflector, for 770 MHz.

Figure 13 .
Figure 13.Signals throughout the processing chain previous to the array signal processing stage.

Figure 15 .
Figure 15.Received power levels provided by the Torrespaña Illuminator of Opportunity (IoO) in the considered radar scenario.

Figure 16 .
Figure 16.Radar scenario located at the roof of the Superior Polytechnic School of Alcalá University.

6. 2 . 1 .
First Array Signal Processing Stage: Target DetectionTaking into consideration the nulls of the main beam associated with the array radiation pattern (Figure12a) and the azimuth coverage area limited by the Televés 4G NOVA beamwidth, the steering angles for the first stage beamformer are: Φ = [−18.66• , 0, 18.66 • ].

Figure 17
presents the NULA azimuth radiation pattern of the multiple simultaneous orthogonal beams.Results show that Alcalá-Meco road and R2 highway are covered for the different orthogonal beams allowing the tracking of cars' trajectories and traffic monitoring, guaranteeing the SLL requirement for all of the steering angles.
(a) Linear representation (b) Polar representation in the considered radar scenario

Figure 18 .
Figure 18.IDEPAR range coverage in the considered radar scenario.(a) Only one Televés 4G NOVA antenna in the surveillance channel; (b) Televés 4G NOVA antenna array in the surveillance channel.

Figure 19 .
Figure 19.CFAR detection performance using a single radiating element (a) and the proposed first array signal processing stage applied to the signals acquired by the designed NULA (b).

Figure 20 .
Figure 20.Targets trajectories in the coverage sector (white).The Polytechnic School of University of Alcalá (EPS, Escuela Politécnica Superior) and a big building made of aluminum (IMMPA, Instituto de Medicina Molecular Príncipe de Asturias) are marked in purple.

Figure 21 .
Figure 21.Comparison of azimuth radiation patterns of NULA with weights calculated for maximizing the directivity (blue) and for guaranteeing SLL < 15 dB (red).

Figure 22 .
Figure 22.DoA based on the beamformer output spectrum for detections in the 30th CPI associated with the considered trajectories.

Figure 23 .
Figure 23.DoA based on the beamformer output spectrum for targets T1, T2 and T3.

Figure 24 .
Figure 24.2D estimated trajectories and GPS data in the radar scenario.

Table 1 .
Ambiguity peaks location in the range domain for the 2 k and 8 k modes, due to Guard Interval (GI), Pilot Carriers (PC) and Transport Parameter Signaling (TPS).

Table 4 .
Main characteristics of the azimuth plane radiation patterns generated with maximum directivity weighting vectors of the ULAs with inter-element spacings of 210 mm and 270 mm and the NULA designed in Section 3.2, for N = 5 commercial antennas.The considered steering angles are φ = [−30, −15, 0, 15, 30] • .

Table 5 .
Main characteristics of the azimuth plane radiation patterns generated with the maximum directivity weighting vector of the ULAs with inter-element spacings of 210 mm and 270 mm and the NULA designed in Section 3.2, for N = 11 commercial antennas.The considered steering angles are φ = [−30, −15, 0, 15, 30] • .

Table 6 .
Main characteristics of the azimuth plane radiation patterns generated with the SLL control weighting vector of the ULAs with inter-element spacings of 210 mm and 270 mm and the NULA designed in Section 3.2, for N = 5 commercial antennas.The considered steering angles are φ = [−30, −15, 0, 15, 30] • .

Table 7 .
Main characteristics of the azimuth plane radiation patterns generated with the SLL control weighting vector of the ULAs with inter-element spacings of 210 mm and 270 mm and the NULA designed in Section 3.2, for N = 11 commercial antennas.The considered steering angles are φ = [−30, −15, 0, 15, 30] • .

Table 9 .
Representative average values of the estimated BRCS of the car model selected for the study.BRCS, Bistatic Radar Cross-Section.

Table 10 .
Detection performances evaluated for the controlled car.

Table 11 .
Number of detections with respect to the number of Coherent Processing Intervals (CPIs) associated with the considered trajectories.

Table 12 .
Signal-to-Interference Ratios (SIR) study for the selected targets.