Passive radar systems can detect targets by utilizing readily available, non-cooperative illuminators of opportunity (IOs) [1, 2, 3, 4], and possess a number of benefits compared with active radar systems. In particular, a passive radar is smaller and less expensive because it does not need a transmitter. And many IOs (e.g., cellular base stations , digital television broadcasting , digital audio broadcasting (DAB) ) are available for passive sensing, as such a passive radar can operate without causing interference to existing communication systems.
A main challenge associated with passive radar is that the IOs are non-cooperative, and the transmitted signals are unknown and not under control. Hence, the conventional matched filter cannot be easily implemented. In addition, the direct-path signal is much stronger than the target reflections, making it difficult to detect and track targets. To solve those problems, a passive radar usually makes use of an additional separate channel, referred to as the reference channel, to collect the transmitted signal in order to eliminate the direct-path signal and clutters in the surveillance channels (SCs) [8, 9]. In addition, the reference signal can also be used to implement approximate matched filtering to detect targets. However, the reference signal is noisy and the target detection performance is usually significantly degraded [1, 3].
Orthogonal frequency-division multiplexing (OFDM) techniques are widely employed in many modern wireless communication systems, e.g., 4G wireless cellular , digital video broadcasting (DVB) , DAB  and wireless local area network (LAN) [12, 13]. For an OFDM passive radar system, demodulation can be implemented by using the reference signal [14, 15, 3]. Since demodulation provides better accuracy than directly using the reference signal, a more accurate matched filter can be implemented based on the estimated data symbols, and the performance of passive radar can be greatly improved. Moreover, the reference channel is not always necessary because data symbols can also be directly demodulated based on the received signal in SC.
Some target detection algorithms using OFDM signals have been proposed [3, 11, 15, 16, 12]. In , a method for detecting a moving target in the presence of multi-path reflections is proposed based on adaptive OFDM radar. In [15, 11, 12], by assuming that the demodulation is perfect, the delays and Doppler shifts of targets are estimated based on matched filtering. And in 
, the MUltiple SIgnal Classifier (MUSIC) and the compressed sensing (CS) techniques are employed to obtain a better target resolution and clutter removal performance. In , by using the received OFDM signal from an un-coordinated but synchronizated illuminator, a delay and Doppler shift estimation algorithm is proposed taking into account the demodulation error. The atomic norm (AN) is used to enforce the signal sparsity in the delay-Doppler plane and the -norm is used to enforce the sparsity of the demodulation error signal. Then, a convex semidefinite program (SDP) is solved to obtain the estimate of the target delays and Doppler shifts.
The present contribution is aimed at extending the results of  by using multiple receivers to achieve the target position and velocity estimation based on the OFDM signal emitted by a totally un-coordinated and un-synchronizated illuminator. Assuming that data demodulation is performed separately by a communication receiver based on the direct-path signal, and the error-prone estimated data symbols are made available to the passive radar receivers. Then, a two-stage procedure is proposed to estimate the positions and velocities of targets.
The first stage is aimed at estimating the delay-Doppler shift and demodulation error by exploiting two types of sparsity: on one hand, as targets and clutters are sparsely distributed in space, the reflected signals hitting the radar receivers are sparse; on the other, the demodulation error rate of a communication system is typically low under normal operating conditions and hence the demodulation error signal is also sparse. Since the delays and Doppler shifts of the targets are continuous parameters, conventional CS tools  may lead to unsatisfactory performance  when the signals cannot be sparsely represented by a finite discrete dictionary [20, 21, 22]. We make use of the recently developed mathematical theory of continuous sparse recovery for super-resolution [23, 24, 25], and especially the AN minimization techniques which have been successfully applied for continuous frequency recovery, line spectral estimation and direction-of-arrival estimation [25, 26, 27]. Note that, unlike the convex problem of the delay-Doppler estimation in  for one receiver, in our model, different receivers share the same estimated data symbols and impose the same constraint, which yields a non-convex problem due to existence of the produce term of decision variables. Hence, we use non-convex factorization (NF) to transform the problem to a smooth unconstrained optimization problem, which is then solved by a conjugate gradient descent (CGD) algorithm.
The second stage is aimed at determining the target positions and velocities based on the estimates in the first stage. Since the illuminator and receivers are un-synchronizated, we utilize the delay differences between different receivers to calculate each target position. The first method numerically solves a set of nonlinear equations, and the second method utilizes the back propagation (BP) neural network [28, 29, 30] to estimate the target position, which is more computationally efficient. The corresponding target velocity can then be determined based on the estimated target position and Doppler shift. Extensive simulation results are provided to illustrate that the proposed methodology can estimate the target positions and velocities accurately.
The remainder of the paper is organized as follows. In Section II, we present the signal model of the OFDM passive radar and set up the problem. In Section III, we develop a delay-Doppler estimator based on conjugate gradient descent. In Section IV, we discuss methods for estimating the locations and velocities. Simulation results are presented in Section V. Section VI concludes the paper.
Ii System Descriptions & Problem Formulation
Ii-a System Descriptions
As shown in Fig. 1, we consider a passive radar system consisting of () receivers and one non-cooperative illuminator, that aims to estimate the locations and velocities of multiple targets in a three-dimensional cartesian coordinate system. Suppose that the coordinates of the illuminator and receiver are and , respectively. Assume that there are reflectors in the surveillance area, include targets and clutters. Note that we consider clutters as zero-velocity targets. Let and be the location and velocity of the -th reflector, respectively. Then, the traveling time from the illuminator to the -th receiving antenna due to the -th reflector is
where is the speed of light in free-space; denotes the -norm. And the corresponding Doppler shift is given by 
where denotes the wavelength of the carrier.
Due to the reflections of targets and clutters, the received signal at the -th receiver is given by
where denotes the convolution operator; is the unit impulse function; is the -th path’s complex gain at the -th receiving antenna; is the unknown communication signal; is the direct-path (illuminator-to-receiver) signal111Note that the direct-path signal can be suppressed by the spatial filtering method in  or by using a reference channel to collect the direct-path signal [1, 2, 32], i.e., using a narrow beam antenna towards the transmitter to receive the direct-path signal , or using the side-lobe of the receiving antenna to receive the direct-path signal .; is the synchronization error between the transmitter and receivers; and is a white, complex circularly symmetric Gaussian process.
Ii-B OFDM-based Passive Radar Signal Model
In this paper, we assume that the signal is the OFDM signal that is widely adopted in contemporary wireless communication systems. The OFDM system consists of sub-carriers, with data sub-carriers and cyclic prefix (CP) sub-carriers. The duration of an OFDM block is with being the “sub-pulse duration.” Then, the transmitted baseband OFDM signal over blocks is given by
where is the -th normalized data symbol block, such that with denoting the complex conjugate operator; and
At each radar receiver , suppose that the direct-path signal is first removed, and we only refer to the baseband signals by assuming that down-conversion has been performed. The CP is removed assuming that its length is no less than the maximum path delay, i.e., . Note that the data symbols are unknown, but can be estimated by demodulation using the direct-path signals . However, demodulation may be error-prone. Hence in the following we assume that an estimate of the data symbols, , is available such that
where denotes the corresponding demodulation error. Furthermore, we assume that the velocity of the target is low, such that . Hence the phase rotation due to the Doppler shift can be approximated as constant over an OFDM block duration , i.e., [15, 3]
At each receiver , in the -th OFDM block matched filtering is performed to obtain, for ,
Ii-C Problem Formulation
Let us now define
and the steering vectors
Correspondingly, the response matrices are defined as
Then (9) can be written as the following matrix form
where denotes the Hadamard product; denotes the diagonal matrix whose diagonal entries are ; , , and are matrices whose -th element are , , and , respectively.
Further denote , , and
and being the Kronecker product. Then we vectorize in (20) to obtain
where is the Khatri-Rao product; is a matrix whose -th column has the form of . Finally, (23) can be rewritten in the following matrix form
where the -th columns of , and are , and , respectively.
In this paper, we first estimate the delays and Doppler shifts contained in from the received signals . Then based on these estimates, we further estimate the locations and velocities of the reflectors , and those with are considered clutters.
Iii Stage 1: Delay-Doppler Estimation
In this section, we propose a CGD method to estimate the delays and Doppler shifts in (24). We first formulate a non-convex optimization problem by exploiting two types of sparsity. Then we relax the non-convex optimization problem to a smooth unconstrained form. The smoothed problem can then be solved via CGD.
Iii-a Non-convex Optimization Problem Setup
We will exploit the following two types of sparsity: firstly, the number of reflectors in (21); secondly, assuming that the demodulation error rate is low, then has a small number of non-zero entries, i.e., with being the -norm. Since the delays and the Doppler shifts take continuous values, the atomic norm [25, 33] is used to exploit the first type of sparsity. Let be the set of atoms, where is defined in (22). Then the 2D atomic norm  associated to for is defined as
The atomic norm can enforce sparsity in the atom set . Note that columns in are independent with their own sparsities. On this basis, our delay-Doppler estimation problem can be formulated according to (24) as:
where denotes the -norm, and are the weight factors.
However, finding the harmonic components via atomic norm is an infinite programming problem over all feasible and . For the convenience of calculation, we use the following equivalent form of (25) for [25, 34]
where denotes the trace operator, stands for a positive semidefinite matrix, and takes as input a matrix
and outputs an block Toeplitz matrix
where denotes the Toeplitz matrix whose first column is the last elements of the input vector. More specifically, we have
Iii-B Conjugate Gradient Descent Algorithm
Then problem (34) is rewritten as
where denotes an inverse operation on the input block Toeplitz matrix, and outputs a matrix. In particular, if we partition the block Toeplitz matrix into blocks, i.e.,
then the -th element of is given by
To solve (36) via the CGD algorithm, we need relax it to a smooth unconstrained form. Hence, we first replace the constraint by the penalty term , and approximate the -norm by a twice continuously differentiable function 
where denotes the -th element in and is a weight parameter, which controls the smoothing level.
Furthermore, we remove the constraints by setting with , such that is chosen minimally according to the ranks of . In particular, we have the following lemma. The proof is given in Appendix A.
The above Lemma 1 shows that each in the solution to (36) should be rank-. It is mentioned in  that an algorithm can be accelerated through non-convex factorization if the solution is low-rank, since the size of the optimization variables is significantly reduced (see Fig. 2). Hence, let be the upper bound on the number of reflectors. We can then let such that the constraints and are both satisfied. Then, (36) can be rewritten as the following smooth unconstrained optimization problem
The CGD algorithm for solving (42) performs the following iterations
where is the step size, which is chosen according to the backtracking line search , given in Appendix B, to guarantee that the objective function does not increase with ; and are the search directions
with and , where and are the gradients, which are derived in Appendix C; and