2-D Coherent Integration Processing and Detecting of Aircrafts Using GNSS-Based Passive Radar

Abstract

Long time coherent integration is a vital method for improving the detection ability of global navigation satellite system (GNSS)-based passive radar, because the GNSS signal is not radar-designed and its power level is very low. For aircraft detection, the large range cell migration (RCM) and Doppler frequency migration (DFM) will seriously affect the coherent processing of azimuth signals, and the traditional range match filter will also be mismatched due to the Doppler-intolerant characteristic of GNSS signals. Accordingly, the energy loss of 2-dimensional (2-D) coherent processing is inevitable in traditional methods. In this paper, a novel 2-D coherent integration processing and algorithm for aircraft target detection is proposed. For azimuth processing, a modified Radon Fourier Transform (RFT) with range-walk removal and Doppler rate estimation is performed. In respect to range compression, a modified matched filter with a shifting Doppler is applied. As a result, the signal will be accurately focused in the range-Doppler domain, and a sufficiently high SNR can be obtained for aircraft detection with a moving target detector. Numerical simulations demonstrate that the range-Doppler parameters of an aircraft target can be obtained, and the position and velocity of the aircraft can be estimated accurately by multiple observation geometries due to abundant GNSS resources. The experimental results also illustrate that the blind Doppler sidelobe is suppressed effectively and the proposed algorithm has a good performance even in the presence of Doppler ambiguity.

1. Introduction

Over the past few years, passive radar has attracted more and more attention and developed very quickly, due to its low cost of operation and maintenance, no need for frequency allocations, and difficulty of jamming [1,2]. Much research about FM, WIFI, satellite TV, and global navigation satellite system (GNSS)-based passive radar have been carried out theoretically and practically [3,4,5,6,7,8]. The common problem of passive radar is that the power level is not high enough for detecting target effectively, and the GNSS signal is no exception, because those transmitted signals are not designed for radar application [9]. However, compared with other illuminators, GNSS signal is widely used recently because of its global coverage and plentiful satellite resources.

To overcome the challenge associated with low power level of the transmitted signals in GNSS-based passive radar, the target forward scattering effect can be taken into consideration [10,11,12]. When the target is crossing the baseline, its radar cross section (RCS) will dramatically increase, and therefore it will be more easily detected in this case. However, the unique geometry will limit its wide application, and sometimes, there is even no available GNSS satellite for forward scattering. For normal bistatic geometry in GNSS-based passive radar, a high receiver antenna gain with extremely long integration duration is highly recommended [13], and the target detection applications are normally focused on medium- or large-sized targets, like aircrafts or marine targets [6,13]. Then, the accompanying problem is long duration, 2-dimensional (2-D), coherent integration processing. For static targets in GNSS-based passive radar, 2-D coherent gain can be easily obtained by synthetic aperture radar (SAR) imaging algorithms [14,15,16]. However, for air moving targets, like an aircraft, it is really difficult to perform coherent processing in either the azimuth or the range dimension with unknown target motion parameters. Due to large Doppler frequency migration (DFM) and across range unit (ARU) effect caused by the fast motion of aircraft, it is not easy to implement coherent processing in the azimuth direction [17,18]. Besides, since GNSS signal is Doppler-intolerant [19], the GNSS-based passive radar is a kind of phase-coded radar and a range-matched filter with target Doppler frequency is necessary to obtain the range coherent gain (see Appendix A). In traditional passive radar, the cross-ambiguity function (cross-correlation) is widely used to obtain the target’s range Doppler parameter or delay Doppler map (DDM) [20,21,22]. However, the cross-ambiguity function method is only applicable in short coherent time case or the Doppler spreading phenomenon could be ignored.

For azimuth coherent processing in the moving target case, the moving target detection (MTD) method was firstly proposed when the DFM and ARU effect can be ignored [23]. Then, the keystone transform (KT) was proposed for range cell migration compensation (RCMC) by reformatting the echo signal [24,25,26], but the computational complexity of this method is very high. For low-complexity blind compensation, a Hough transform-based integration method was proposed by transforming the signal into the range-Doppler-time space [27]. However, it cannot be applied when the SNR is extremely low. Recently, the Radon Fourier Transform (RFT) was proposed to carry out long time coherent processing via joint searching along the range and velocity dimensions [18,28,29,30], and the blind speed effect can also be suppressed by weighted window processing [31]. However, integration loss may exist because of quadratic phase error. Based on the principle of RFT, a generalized and complex method called Radon-fraction Fourier transform was proposed for maneuvering target detection [32]. All those research works demonstrate that the RFT operation is a useful tool for moving target long time coherent integration. For Range-matched filtering in phase-coded radar, the traditional method is searching for the peak with different Doppler shift, which is very time consuming [33]. Besides, it may not work, owing to the low power level of echo signal in GNSS-based passive radar.

In this paper, motivated by the previous work on RFT, a novel 2-D coherent integration processing and detecting algorithm is proposed to enhance the performance of GNSS-based passive radar for aircraft detection. Firstly, the phase term and range-walk caused by the motion of GNSS is compensated. Then, based on the principle of RFT, azimuth coherent integration processing is performed with range-walk removal and residual quadratic phase error (QPE) estimation and compensation. In addition, the fast implementation method for azimuth coherent integration processing is presented based on the Chirp-Z transform. After that, a modified Doppler shifting matched filter is proposed for range compression. This is followed by aircraft detection by the traditional track and detection (TAD) [34] or track before detection (TBD) method [35], with the focused range-Doppler data after 2-D coherent integration processing.

The paper is organized as follows. Analysis of Doppler characteristics of air targets and its echo model are presented in Section 2. In Section 3, the proposed 2-D coherent integration processing and target detection algorithm are derived in detail, and its processing framework is also presented. In Section 4, numerical simulations are carried out to demonstrate the validity and feasibility of our proposed algorithm, and the aircrafts’ motion parameters, including position and velocity, can be obtained accurately by multi-static joint experimental results; besides, the performance of the proposed algorithm is analyzed and discussed, from three aspects: System power budget and available coherent time for aircraft detection, best choice of the Doppler rate searching step, and processing performance when Doppler ambiguity exists. The last section concludes the paper and presents its future research directions.

2. Air Target Echo Model in GNSS-Based Passive Radar

First of all, the geometry of GNSS-based passive radar will be described and the characteristics of system range and Doppler history analyzed based on the acquisition geometry, followed by derivation of the air target echo model with GNSS signals.

2.1. GNSS-Based Passive Radar Geometry and Doppler Characteristic Analysis

The air target detection geometric configuration of GNSS-based passive radar is shown in Figure 1, where the radar receiver is at the origin of the Cartesian coordinate system, while the X and Y axes point to the East and North directions, respectively. The air target is located at P_tar(X_A,Y_A,Z_A) and moves at a constant velocity V_tar(V_x, V_y, 0) and in a constant level. The GNSS transmitter is located at P_Tran(X_G,Y_G,Z_G) with a velocity vector V_Tran(V_Tx, V_Ty, V_Tz).

Based on the system geometry, the range history is the sum of the instantaneous receiver slant range R_R(η) and transmitter slant range R_T(η), which can be rewritten by Taylor series expansion.

$$R\left(\eta \right)=R_R\left(\eta \right)+R_T\left(\eta \right)=R_{ref}+\lambda f_d\eta +\frac{\lambda }{2}{f}_r{\eta }^2$$

(1)

where η, R_ref, λ, f_d, and f_r denote the illumination time, reference range, signal wavelength, target Doppler centroid, and Doppler modulated rate, and f_d and f_r are derived as follows [36]:

$$f_d=\frac{1}{\lambda }{\frac{dR\left(\eta \right)}{d\eta }|}_{\eta =0}=f_{d0}+f_{d,v},f_r=\frac{1}{\lambda }{\frac{d^{2}R\left(\eta \right)}{d{\eta }^2}|}_{\eta =0}=f_{r0}+f_{r,v}$$

(2)

where η = 0 means the illumination centroid time, f_d0 and f_r0 represent the Doppler centroid frequency and Doppler frequency modulated rate of the stationary target, and f_d,v and f_r,v are the corresponding variations caused by target motion. Besides, in GNSS-based passive radar, the reference signal, and its Doppler are also important for signal synchronization and range compression. As both the transmitter and the receiver co-ordinates are known, it is easy to remove the Doppler of reference signal from the synchronization outputs [14]. Therefore, to simplify the subsequent derivation, the Doppler of reference signal are not considered in this manuscript.

The additional unknown f_d,v and f_r,v will cause Doppler spectrum shift and image defocusing, which is a serious problem for long duration coherent integration processing and target detection. Assuming the air target is moving in the XOY plane with a velocity 220 m/s, Figure 2 shows the variations of f_d,v and f_r,v, using the parameters listed in

Table 1. Main simulation parameters.

Parameters	Values	Parameters	Values
Wavelength	0.255 m	SVN #2 position	(1.0234, −1.5530,1.2411) × 104 km
Sampling rate	42.0 MHz	SVN #7 position	(0.9767, 1.2291, 1.4880) × 104 km
$\Delta f_v$	0.2 Hz	SVN #12 position	(−1.8017, 1.6150, 0.4690) × 104 km
$\Delta {\widehat{f}}_r$	0.04 Hz/s	SVN #2 velocity	(187.9, −2113.8, −1799.3) m/s
NV	15,000	SVN #7 velocity	(−2398.5, −632.5, 1426.5) m/s
fp	1000 Hz	SVN #12 velocity	(−2097.3, −729.0, −2321.4) m/s
Experiment scenario Center	(60, 40, 10) km	Target RCS [41]	20 dB
Antenna gain (Receiver)	35 dB	Coherent time	5 s
Target 1 position	(70, 20, 10) km	Target 1 velocity	(150, 120) m/s
Target 2 position	(65, 30, 10) km	Target 2 velocity	(−135, 140) m/s
Target 3 position	(50, 50, 10) km	Target 3 velocity	(100, −170) m/s

. As shown in Figure 2, both f_d,v and f_r,v are greatly influenced by passive radar geometry and target motion direction. f_d,v is mainly determined by the equivalent range velocity, and f_r,v is mainly influenced by the equivalent azimuth velocity [36]. Normally, as depicted in Figure 2a, the range of f_d,v is from −1500 Hz to 1500 Hz when the target velocity is smaller than 220 m/s (approximately 800 km/h). However, the pulse repetition frequency (PRF), f_p, of the GNSS signal is 1000 Hz. As a result, the Doppler ambiguity phenomenon should be considered in the signal processing; otherwise, the blind speed will be inevitably detected. On the other side, f_r,v is about 2.6 ± 0.2 Hz/s on this condition, which is relatively small because of the air target is far away from the receiver, and the corresponding range-curvature term is smaller than one range cell normally.

2.2. Air Target Echo Signal Model of GNSS-Based Passive Radar

In a GNSS system, the transmitted signal is a code-modulated continuous wave, the echo signal of GNSS-based passive radar forms a one-dimensional vector, and the target reflected GNSS signals from all visible satellites are received by the receiver antenna. According to [16], after quadrature demodulation, the target echo signal can be modeled as:

$$S\left(\eta \right)={\displaystyle \sum _{i=1}^I{\displaystyle \sum _{n=-\infty }^{\infty }{\sigma }^i{A}_{T}Q\left[\eta -n{T}_p-\frac{R^i\left(\eta \right)}{c}\right]}\exp \left\{-j\frac{2\pi R^i\left(\eta \right)}{\lambda }\right\}}$$

(3)

where superscript i represents using the i-th visible satellite as the transmitter, σ is the target radar cross-section, A_T is the signal amplitude, Q(·) is the modulation code, and T_p, c, and λ denote signal period, speed of light, and signal wavelength, respectively. Then, set η = nT_p + τ, t_n = nT_p, and the range history can be written as:

$$R^i\left(\eta \right)=R^i\left(n{T}_p+\tau \right)\approx R^i\left(t_n\right)+\lambda f_d^i\tau$$

(4)

The Doppler frequency is also considered during the inter-period time, which is totally different with the traditional pulse radar. Besides, the second-term with inter-period time is ignored in (4).

Therefore, based on (4), the echo signal of GNSS-based passive radar can be rewritten in a 2-D form.

$$S\left(t,\tau \right)={\displaystyle \sum _{i=1}^I{\sigma }^i{A}_{T}Q\left[\tau -\frac{R^i\left(t\right)+\lambda f_d^i\tau }{c}\right]\exp \left\{-j\frac{2\pi R^i\left(t\right)}{\lambda }\right\}}\cdot \exp \left\{-j2\pi f_d^i\tau \right\}$$

(5)

where the subscript n in t_n is ignored for simplicity. As the transmitted code in each GNSS satellite is high-rate pseudo random noise (PRN) sequence, it has particularly excellent auto-correlation and cross-correlation properties, and these cross-correlation values are so small that they usually can be ignored. As a result, corresponding PRN code could not only be used as the matched filter but also be employed for separating the echo signal from different visible GNSS satellites. Hence, to simplify the subsequent derivation, one specific visible GNSS satellite is chosen as the transmitter, and the echo signal can be written as:

$$S\left(t,\tau \right)=\sigma A_{T}Q\left[\tau -\frac{\left(R_{ref}+\lambda f_d^{}t\right)}{c}\right]\exp \left\{-j\frac{2\pi R\left(t\right)}{\lambda }\right\}\exp \left\{-j2\pi f_d^{}\tau \right\}$$

(6)

where R_ref is the target reference range in center time, and the range-curvature in the envelop term is ignored, as it is much smaller than one range cell.

3. The Proposed 2-D Coherent Integration Processing and Target Detection Algorithm

Compared with the traditional radar signal, the signal power of GNSS is extremely weak. Consequently, to improve the performance of GNSS-based passive radar, not only the range integration gain, but also azimuth integration gain with long duration, must be obtained. In this section, a novel 2-D coherent integration algorithm is proposed, and based on this, a framework for aircraft target detection is presented. Firstly, the principle of RFT in introduced in Section 3.1. Then, the proposed algorithm is fully detailed in Section 3.2, with its fast implementation method and the target detection framework.

3.1. The Principle of RFT

For moving targets in radar applications, the echo signal after range compression s(t, r) is distributed on the range history line R(t) = v₀t + r₀, where r₀ is the initial target position and v₀ is target radial velocity. To coherently integrate the moving target’s azimuth echoes pulse by pulse, the traditional RFT for radar signal is by jointly searching range and velocity (r, v), which is [18]:

$$R\left(v,r\right)={\displaystyle {\int }_{-\infty }^{\infty }s\left(t,\frac{2\left(vt+r\right)}{c}\right)\exp \left\{j\frac{4\pi vt}{\lambda }\right\}}dt$$

(7)

As Doppler frequency is f_v = 2v/λ, the range history can be rewritten by R(t) = λf_v0t/2 + r₀, where f_v0 is the target Doppler frequency. So, the RFT in the range-Doppler plane is:

$$R\left(f_v,r\right)={\displaystyle {\int }_{-\infty }^{\infty }s\left(t,\frac{2r+\lambda f_{v}t}{c}\right)\exp \left\{j2\pi f_{v}t\right\}}dt$$

(8)

Then, by defining the discrete searching range and Doppler parameters, the discrete form of traditional RFT has also been defined as [18]:

$$\begin{array}{c}R\left(i,n\right)={\displaystyle \sum _{m=1}^{N_a}s\left[m{T}_p,\frac{2r\left(n\right)+\lambda \cdot f_v\left(i\right)\cdot m{T}_p}{c}\right]\exp \left\{j2\pi \cdot f_v\left(i\right)\cdot m{T}_p\right\}}\\ r\left(n\right)=r_{\min }+n\Delta r,n=1,2,3,\dots ,N_r,f_v\left(i\right)=f_{\min }+i\Delta f_v,i=1,2,3,\dots ,N_f\end{array}$$

(9)

where Δr and Δf_v are range and Doppler searching steps, r_min and f_min are the minimum target range and Doppler searching value, and N_a, N_r, and N_f are the sample number of azimuth, range, and Doppler domains, respectively. As shown in (7)–(9), the RFT actually performs echo signal integration along a parameterized line defined by (r, v) or (r, f_v). When the parameterized line (r₀, v₀) or (r₀, f_v0), coincides with the target range history, all echoes distributed on the line are integrated to a peak, and full coherent-integration gain processing is achieved.

In the traditional RFT, the Doppler rate or target high-order motion isn’t considered, and the monostatic radar is adopted. Then, the generalized RFT [28] has also been proposed using curvilinear integral for high-order motion detection, but higher dimensional parameter searching is needed, which is time and memory consuming.

3.2. Air Target Detection Based on 2-D Coherent Integration

To guarantee the effectiveness and accuracy of air target detection, a sufficiently high SNR is crucial, which means long time illumination and 2-D coherent integration are needed. In the proposed novel air target detection algorithm, range-matched filtering and Doppler searching are combined to achieve the 2-D coherent-integration gain. Based on the deduced novel 2-D coherent integration algorithm, the framework of target detection using GNSS-based passive radar is shown in Figure 3, which mainly has four parts: The first is signal preprocessing to compensate the phase term and range-walk caused by the motion of GNSS satellite; the second is azimuth coherent integration processing with QPE compensation and range-walk removal; the third is range compression with a modified Doppler shift matched filter, after which the echo signal is totally focused on the range-Doppler plane; the last is air target detection based on the TAD or TBD method. The following provides details of the operations.

3.2.1. Signal Preprocessing

The signal preprocessing stage is performed to compensate the phase term and remove the range-walk, and both of them are caused by GNSS satellite’s motion. It starts with a range FFT to transform the echo signal into the range frequency domain, which is:

$$\begin{array}{ll}{S}_r\left(t,f_{\tau }\right)& ={\displaystyle \int S\left(t,\tau \right)\exp \left\{-j2\pi f_{\tau }\tau \right\}d\tau }\\ & \approx A{Q}_f\left(f_{\tau }+f_d^{}\right)\exp \left\{-j\frac{2\pi f_{\tau }{R}_{ref}}{c}\right\}\cdot \exp \left\{-j\frac{2\pi f_{\tau }\lambda f_{d}t}{c}\right\}\cdot \exp \left\{-j2\pi f_{d}t\right\}\exp \left\{-j\pi f_r{t}^2\right\}\end{array}$$

(10)

where Q_f(·) is the FFT form of the GNSS modulation code,$A=\sigma A_T\exp \left\{j\frac{2\pi R_{ref}}{\lambda }\right\}$. Then, the phase term and range-walk caused by f_d0 and f_r0 are compensated with the 2-D filter $H_p\left(t,f_{\tau }\right)$.

$$H_p\left(t,f_{\tau }\right)=\exp \left\{j\pi \left(f_{r0}{t}^2+2f_{d0}t\right)\right\}\cdot \exp \left\{j\frac{2\pi f_{\tau }\lambda f_{d0}t}{c}\right\}$$

(11)

In (11), the f_d0 and f_r0 are unknown, because the target position is unknown in advance. However, as the satellites orbit is about 20,000 km, the changing of slant range R_T in different range cell is very small compared with the R_T value itself. So, the difference of f_d0 and f_r0 in each position of the whole experiment scenario is quite small, and any range cell’s f_d0 and f_r0 is enough to guarantee the accuracy of 2-D filter in (11). In this situation, it is recommended to use the f_d0 and f_r0 in the scene center to construct the 2-D filter.

Then, the signal becomes:

$$S_p\left(t,f_{\tau }\right)=S_r\left(t,f_{\tau }\right)\cdot H_p\left(t,f_{\tau }\right)=A{Q}_f\left(f_{\tau }+f_d^{}\right)\exp \left\{-j\frac{2\pi f_{\tau }{R}_{ref}}{c}\right\}\cdot \exp \left\{-j\frac{2\pi f_{\tau }\lambda f_{d,v}t}{c}\right\}\exp \left\{-j\pi \left(2f_{d,v}t+f_{r,v}{t}^2\right)\right\}$$

(12)

3.2.2. Azimuth Coherent Integration

To obtain the coherent integration processing gain in azimuth, the principle of RFT is adopted at this stage. After the signal preprocessing, there is still residual quadratic phase error (QPE) caused by the unknown, f_r,v, and residual range-walk term caused by the unknown, f_d,v. For QPE estimation and compensation, many methods have been proposed in SAR imaging, and the simple and easy way is based on amplitude auto-focus method [37]. Besides, the residual range-walk can be removed by Doppler searching. Therefore, the azimuth coherent integration processing can be performed by QPE compensation, range-walk removal, and RFT, which is:

$$\begin{array}{ll}{S}_{2D}\left(f_v,f_{\tau }\right)& ={\displaystyle \int S_p\left(t,f_{\tau }\right)\cdot \exp \left\{j\pi {\widehat{f}}_r{t}^2\right\}\cdot \exp \left\{j\frac{2\pi f_{\tau }\lambda f_{v}t}{c}\right\}\cdot \exp \left\{j2\pi f_{v}t\right\}}dt\\ & =A{Q}_f\left(f_{\tau }+f_d^{}\right)\exp \left\{-j\frac{2\pi f_{\tau }{R}_{ref}}{c}\right\}{\displaystyle \int \exp \left\{j\frac{2\pi f_{\tau }\lambda \left(f_v-f_{d,v}\right)t}{c}\right\}\exp \left\{j2\pi \left(f_v-f_{d,v}\right)t\right\}\exp \left\{j\pi \Delta f_r{t}^2\right\}}dt\end{array}$$

(13)

where

• The first phase term in the first line of Equation (13) is QPE compensation, and ${\widehat{f}}_r$ is the estimated Doppler rate parameter. As shown in Figure 2b, the variation range of Doppler rate is very small, and the maximum variation value is smaller than 0.6 Hz/s. Hence, we can use parallelized step search to estimate the Doppler rate (as shown in Figure 4), and this method has a high searching efficiency because of the small variation range of f_r,v. Normally, the estimated Doppler rate can be expressed as:

  $${\widehat{f}}_r^k={\widehat{f}}_r^0+k\Delta {\widehat{f}}_r,k=0,1,2,\dots ,N_k$$

  (14)

  where ${\widehat{f}}_r^0$, $\Delta {\widehat{f}}_r$ and $N_k$ denote the initial Doppler rate, Doppler rate searching step, and Doppler rate searching number, respectively. Finally, when ${\widehat{f}}_r^k$ matches the true value of target Doppler rate or $\Delta f_r$ equals zero, the azimuth coherent gain could be achieved and the maximum peak value of target will be obtained.
• The second phase term in the first line of Equation (13) is residual range-walk removal, and f_v is the Doppler searching parameter. As shown in Equation (9), non-integer range cell is inevitable during the range history line searching, so range interpolation processing is necessary for traditional RFT. However, in Equation (13), when f_v matches the target’s Doppler frequency, f_d,v, all range-walk term is removed, and the echo signal will distribute along the line perpendicular to the range direction. Therefore, the non-interpolation processing is needed in the proposed algorithm, and the following azimuth integral processing is greatly simplified.
• The last phase term in the first line of Equation (13) is RFT processing, or called azimuth integral. As the range-walk term is removed, the echo signal is distributed along the azimuth direction. So, the RFT operation is similar to azimuth FFT operation.

Therefore, as shown in the azimuth integral in the second line of Equation (13), the azimuth coherent integration gain (see Equation (24)) is achieved when the Doppler searching parameter equals f_d,v, the estimated Doppler rate equals f_r,v, and the azimuth peak will occur at the target Doppler frequency point.

However, in Equation (13), the range-walk removal term is a three-dimensional filter, which means the whole azimuth coherent integration processing will have three nested loops in digital processing, and the computational complexity of this stage is o(N_AN_RN_V), where N_A, N_R, and N_V denote azimuth sample number, range sample number, and Doppler searching number, respectively. Therefore, the azimuth coherent integration processing will take a lot of memory and increase the computational complexity. To improve the algorithm’s efficiency, a fast implementation of azimuth coherent integration based on the chirp-Z transform (CZT) [38] is presented in the following.

Firstly, we define the discrete series of searching Doppler frequency, azimuth time, and range frequency as:

$$\begin{array}{l}{f}_v=i\Delta f_v,i=-\frac{N_V}{2},-\frac{N_V}{2}+1,\dots ,\frac{N_V}{2}-1\\ t=m{T}_p,m=-\frac{N_A}{2},-\frac{N_A}{2}+1,\dots ,\frac{N_A}{2}-1\\ f_{\tau }=n\Delta f_{\tau },n=-\frac{N_R}{2},-\frac{N_R}{2}+1,\dots ,\frac{N_R}{2}-1\end{array}$$

(15)

where $\Delta f_v$ and $\Delta f_{\tau }$ denote Doppler searching step and range frequency cell, respectively. As shown in Figure 2a, the searching Doppler frequency is ambiguous, and the searching Doppler frequency series can be rewritten as:

$$f_v=i\Delta f_v+I\cdot f_p,i=-\frac{N_v}{2},-\frac{N_v}{2}+1,\dots ,\frac{N_v}{2}-1$$

(16)

where I and N_v denote the ambiguous Doppler index and non-ambiguous Doppler searching number. Based on this, Equation (13) can be rewritten as:

$$S_{2D}\left(i,n\right)={\displaystyle \sum _{m=-N_a/2}^{N_A/2-1}{S}_p\left(m,n\right)\cdot \exp \left\{j\pi {\widehat{f}}_r{\left(m{T}_p\right)}^2\right\}\cdot \exp \left\{j2\pi m{T}_p\left(1+\frac{\lambda n\Delta f_{\tau }}{c}\right)\cdot \left(i\Delta f_v+I{f}_p\right)\right\}}={\displaystyle \sum _{m=-N_A/2}^{N_A/2-1}{S}_q\left(m,n\right)\cdot Y^{mi}}$$

(17)

where

$$\begin{array}{c}{S}_q\left(m,n\right)=S_p\left(m,n\right)\cdot \exp \left\{j\pi {\widehat{f}}_r{\left(m{T}_p\right)}^2\right\}\cdot \exp \left\{j2\pi I{f}_{p}m{T}_p\left(1+\frac{\lambda n\Delta f_{\tau }}{c}\right)\right\}\\ Y=\exp \left\{j2\pi \Delta f_v{T}_p\left(1+\frac{n\lambda \Delta f_{\tau }}{c}\right)\right\}\end{array}$$

(18)

Then, using the principle of the CZT, we have:

$$\begin{array}{ll}{S}_{2D}\left(i,n\right)& ={\displaystyle \sum _{m=-N_a/2}^{N_a/2-1}{S}_q\left(m,n\right)\cdot Y^{mi}}={\displaystyle \sum _{m=-N_a/2}^{N_a/2-1}{S}_q\left(m,n\right)\cdot Y^{\frac{1}{2}\left[m^2+i^2-{\left(m-i\right)}^2\right]}}\\ & =Y^{\frac{i^2}{2}}\cdot {\displaystyle \sum _{m=-N_a/2}^{N_a/2-1}\left[S_q\left(m,n\right)\cdot Y^{\frac{m^2}{2}}\right]\cdot Y^{-\frac{{\left(m-i\right)}^2}{2}}}\\ & =Y^{\frac{i^2}{2}}\cdot \left\{\left[S_q\left(m,n\right)Y^{\frac{m^2}{2}}\right]\otimes Y^{-\frac{m^2}{2}}\right\}\end{array}$$

(19)

where $\otimes$ denotes the convolution operation. Besides, based on the relationship between convolution and FFT, Equation (19) can be performed in digital processing as:

$$S_{2D}\left(i,n\right)=Y^{\frac{i^2}{2}}\cdot IDF{T}_m\left\{DF{T}_m\left[S_q\left(m,n\right)Y^{\frac{m^2}{2}}\right]\cdot DF{T}_m\left(Y^{-\frac{m^2}{2}}\right)\right\}$$

(20)

where subscript m means the operation in the azimuth direction, DFT_m(·) and IDFT_m(·) denote the Discrete Fourier Transform and Inverse Discrete Fourier Transform, respectively.

Actually, the calculation of (17) contains three nested loops (Doppler frequency, azimuth time, and range frequency), which is very memory and time consuming. But, as shown in (20) with CZT, the 3-D azimuth coherent integration processing is reduced to 2-D operation, and the computational complexity is o(N_RN_V), which is significantly lower than the 3-D processing. More details about computation complexity analysis can be seen in Chapter 4.2.4. Based on (20), the processing flow of azimuth coherent integration is shown in Figure 5.

3.2.3. Range Compression

After azimuth coherent integration processing, the target echo is focused in the Doppler domain. To maximize the signal-to-noise ratio and obtain fine range resolution of the sensed air target, range compression processing is performed. In tradition, the radar transmitted signal, like LFM, can be used as the matched filter for range compression. However, in GNSS-based passive radar, the transmitted PRN sequence is Doppler-intolerant, and a moving air target will introduce an unknown and non-negligible Doppler frequency f_d,v. This will introduce a time-varying phase in the range direction and change the response of range compression. The traditional cross-ambiguity function could be used for range compression with Doppler and range 2-D joint-searching in each signal period [20,21,22]. But it will very time consuming because of the extremely long azimuth coherent time. Therefore, a special matched filter including Doppler frequency is necessary.

Normally, the matched filter is generated by the reference signal. In this manuscript, the Doppler of the reference signal is ignored, because it is easy to remove with a known baseline range history [14]. Besides, the phase errors caused by atmospheric propagation and clock slippage are also ignored, because those errors in reference signal and target signal can be modelled as having the same error, and they are easily removed after range compression [14].

The signal after azimuth coherent integration is transformed into the 2-D frequency domain, so the modified matched filter can be modeled easily in the 2-D frequency domain, and a shifting Doppler in each Doppler line is included.

$$H_r\left(f_v,f_{\tau }\right)=Q_f\left[f_{\tau }+\left(f_v+f_{d0}\right)\right]=DF{T}_{\tau }\left(Q\left(\tau \right)\exp \left\{-j2\pi \left(f_v+f_{d0}\right)\tau \right\}\right)$$

(21)

where subscript τ means the FFT operation in the fast time direction. By multiplying the modified matched filter $H_r\left(f_v,f_{\tau }\right)$ and performing range IFFT operation, the signal after 2-D coherent processing is:

$$S_o\left(f_v,\tau \right)={\displaystyle \int S_{2\mathrm{D}}\left(f_v,f_{\tau }\right)\cdot H_r^{*}\left(f_v,f_{\tau }\right)\exp \left\{j2\pi f_{\tau }\tau \right\}d}{f}_{\tau }$$

(22)

where the superscript * represents the complex conjugate of $H_r\left(f_v,f_{\tau }\right)$. Then, let $R=c\tau$, and replacing the azimuth time integral with summation, the final result can be expressed as:

$$\begin{array}{ll}{S}_o\left(f_v,R\right)=& A{Q}^C{\displaystyle \sum _{m=-N_A/2}^{N_A/2-1}D\left[R-R_{ref}-\delta \left(f_v-f_{d,v}\right)+\lambda \left(f_v-f_{d,v}\right)m{T}_p\right]}\\ & \cdot \exp \left\{j2\pi \left(f_v-f_{d,v}\right)m{T}_p\right\}\cdot \exp \left\{j\pi \Delta f_r{\left(m{T}_p\right)}^2\right\}\end{array}$$

(23)

where Q^C is the gain of range compression, normally Q^C equals B_w/f_p, and B_w is the signal bandwidth. D(·) is the normalized envelope of the auto-correlation result of GNSS transmitted signal, and D(0) = 1. $\delta \left(f_v-f_{d,v}\right)$ is the factor caused by mismatched Doppler frequency in range compression, and $\delta \left(0\right)=1$, which means the range compression is matched filtering when $f_v=f_{d,v}$. Here, as shown in Figure 4, we use parallelized search to estimate the Doppler rate. Then, when the error of estimated Doppler rate $\Delta f_r$ is zero, the signal $S_o\left(f_{d,v},R_{ref}\right)$ is:

$${S_o\left(f_{d,v},R_{ref}\right)|}_{\Delta f_r=0}=A{N}_A{Q}^C$$

(24)

where azimuth sample number, N_A, equals coherent time multiplied by the signal repetition rate, and N_A is the azimuth coherent integration. Considering the geometry of GNSS-based passive radar, each moving target’s range Doppler parameters are totally different. Therefore, even in multi-targets case, all targets are focused at their own Doppler parameter position.

Clearly, the maximum value will occur at the position of target motion parameter $\left(f_{d,v},R_{ref},f_{r,v}\right)$. Therefore, by comparing the peak value of the processing result with different estimated Doppler rate,${\widehat{f}}_r$, the coherent integration gain can be achieved after the 2-D coherent integration processing.

3.2.4. Target Detection

After the 2-D coherent integration processing, the air target is focused on the range-Doppler plane. Apart from the target’s signal, the system noise and clutter also exist in the focused range-Doppler data. Therefore, target detection is performed to obtain the air target’s motion parameters. Normally, the likelihood ratio function (LRF) is used for target detection, and many methods based on LRF have been presented, which can be divided into two categories: Tracking and detection (TAD), and track before detection (TBD) [34,35]. In TAD, the constant false alarm detection (CFAR) is typically used, and to obtain the detection probability of 80% with a constant false alarm probability 10⁻⁶ in the additive white Gaussian noise background, the minimum signal-to-noise ratio (SNR) is 12.8 dB [18]. In TBD, the particle filter (PF) is commonly used, and it is capable of detecting a moving target with SNR as low as 7 dB [35]. In this paper, the CFAR method is used in the target detection part.

4. Experiments and Discussion

To demonstrate the performance of the proposed algorithm and the feasibility of target detection using GNSS-based passive radar, numerical simulations including comparison experiments and multi-static experiments are carried out in Section 4.1. Then, based on the results, the performance of GNSS-based passive radar and the proposed algorithm is discussed in Section 4.2, which mainly has three parts: System power budget analysis and available coherent time for aircraft target detection, best choice of the Doppler rate searching step, and processing performance when Doppler ambiguity exists.

4.1. Experiments and Results

In the experiments, the newly modernized civil signal GPS L5 is adopted, as it is at least ten times wider in signal bandwidth and four times stronger in terms of signal power than traditional civil signal GPS L1 [39,40]. Hence, both of signal bandwidth and signal power are beneficial to aircraft detection, and GPS L5 is a superior option in GNSS-based passive radar. The motion parameters of three targets and the associated simulation parameters are listed in Table 1. First, a comparison between the tradition RFT and the proposed algorithm is presented to verify the validity of the proposed algorithm. Then, multi-satellites (three satellites) joint experiments are presented for aircraft detection and motion parameters estimation.

4.1.1. Comparing Results Using the Proposed Algorithm and Traditional RFT

To illustrate the performance and feasibility of the proposed 2-D coherent processing algorithm, the processing results of target-3 with traditional RFT [18] and our proposed algorithm are shown in Figure 6. GPS SVN #7 is used as the transmitter, and noise is not considered in echoes. As shown in Figure 6b,d, the processing result using the traditional RFT is totally defocused in the Doppler dimension, as the residual QPE (Doppler rate term) caused by the motion of target. Besides, energy loss exists after range compression, as shown in Figure 6c, because the range match filter in traditional RFT is mismatched for an unknown non-zero Doppler frequency. For the proposed 2-D coherent processing algorithm, both Doppler and range profile are accurately focused as shown in Figure 6a. Furthermore, as listed in

Table 2. True and estimated range and Doppler frequency for the three targets.

			Target-1	Target-2	Target-3
SVN #2	fd,v (Hz)	True	748.87	375.28	−835.72
		Estimated	749.20	375.40	−835.80
	Variation from reference range (km)	True	−17.79	−9.76	10.13
		Estimated	−17.78	−9.76	10.13
SVN #7	fd,v (Hz)	True	156.66	−322.39	9.49
		Estimated	156.40	−322.20	9.60
	Variation from reference range (m)	True	7.55	2.92	−2.56
		Estimated	7.55	2.92	−2.56
SVN #12	fd,v (Hz)	True	812.52	−995.59	531.07
		Estimated	812.60	−998.60	531.20
	Variation from reference range (m)	True	21.08	9.68	−15.25
		Estimated	21.08	9.68	−15.24

, the peak of them is equal to target’s Doppler and range parameters, and the accuracy is better than one Doppler searching or range bin. Thus, the proposed method has a good performance on moving target focusing.

4.1.2. Multi-Static Experiments and Aircrafts Motion Parameters Estimation

Due to the abundance of GNSS signal resources, multi-angle illumination from over 20 satellites in any given area is attainable. So, multi-satellite (SVN #2, #7, and #12) joint experiments are carried out to verify the applicability and generality of the proposed 2-D coherent processing and target detection method. With different geometries for those three transmitters, each target’s Doppler and range parameters are totally different. Even the same target’s Doppler and range parameters are changed with different transmitters. All motion parameters of those three targets are listed in Table 2. According to the power budget, those three targets’ ideal SNRs are 13.6, 13.8, and 13.9 dB, respectively. By applying the proposed 2-D coherent processing algorithm, those three targets are well-focused as demonstrated in Figure 7, and their measured SNR is 13.5, 13.6 and 13.8 dB respectively. The energy loss is mainly because of the residual estimated Doppler rate. Even so, the energy loss is small enough and acceptable. Then, based on the 2-D focused result in the range-Doppler plane, target detection experiments are carried out with the CFAR method, and the detection results are listed in Table 2. As shown, all targets have been detected and accurately estimated Doppler and range parameters are obtained.

The connection between a target’s location, velocity, and range Doppler parameters is:

$$\{\begin{array}{l}{R}_{ref}=F\left(X_A,Y_A,Z_A\right)\\ f_{d,v}=G\left(X_A,Y_A,Z_A,V_X,V_Y,V_Z\right)\end{array}$$

(25)

where the functions F(·) and G(·) are determined by (1) and (2), respectively. According to (25), it is impossible to achieve a target’s location and velocity with only a group of range Doppler parameters. Fortunately, as shown in Figure 7 and Table 2, a set of different range Doppler parameters could be achieved with each available GNSS satellite. Therefore, multi-element nonlinear equations can be obtained by using several satellites (at least three to obtain the target’s 3-D motion parameters). Finally, as listed in

Table 3. True and estimated location and velocity of the three targets.

		Target-1	Target-2	Target-3
Target Location (XA,YA,ZA) (km)	True	(70, 20, 10)	(65, 30, 10)	(50, 50, 10)
	Estimated	(70.02, 19.99, 10.00)	(64.98, 30.01, 10.01)	(50.01, 50.02, 10.00)
Target Velocity (Vx,Vy,0) (m/s)	True	(150, 120, 0)	(−135, 140, 0)	(100, −170, 0)
	Estimated	(150.2, 119.9, 0)	(−135.0, 140.1, 0.0)	(99.9, −170.1, 0)

, all targets’ locations and velocities are estimated accurately.

4.2. Discussion

4.2.1. Power Budget and Coherent Time

The system thermal noise always exists in the radar echo signal, and, according to the radar function, the SNR_in of echo signal is:

$$SN{R}_{in}=\frac{P_E{G}_R{\sigma }_T}{R_R^2\cdot k{T}_s{B}_w\cdot F_{n}L}$$

(26)

where P_E is the GNSS signal power near the Earth’s surface, and it equals −154 dB in GPS L5 signal case [40]; G_R is the antenna gain of receiver, and it is set to 35 dB as listed in Table 1; σ is target RCS, for a medium-sized aircraft (like Boeing 737–400), its RCS exceeds 20 dB in L-band, and in some azimuth angle, the RCS is larger than 26 dB [41]; R_R is the range between target and receiver, k is the Boltzmann constant, T_s is system Kelvin temperature, B_w is the signal bandwidth, and F_nL is system noise coefficient. For GNSS-based passive radar with the GPS-L5 signal, SNR_in is roughly smaller than −70 dB. Hence, the 2-D coherent integration gain is necessary to obtain a sufficiently high SNR for target detection.

After the 2-D coherent integration processing, the azimuth and range processing gain is achieved, and the SNR_o of processing result is:

$$SN{R}_o=SN{R}_{in}\cdot T_p{B}_w\cdot f_p{T}_c=SN{R}_{in}\cdot B_w\cdot T_c$$

(27)

where T_c is the coherent time, T_pB_w and f_pT_c are the range and azimuth processing gain, respectively. Normally, to guarantee the effectiveness of target detection, the minimum SNR (SNR_min) of the processing results should be obtained. So, the coherent time, T_c, needs to be longer than the lower bound, which is:

$$T_c\ge \frac{SN{R}_{\min }}{SN{R}_{in}{T}_p{B}_w{f}_p}$$

(28)

Furthermore, during the derivation of the proposed algorithm, the assumption that the range-curvature term is less than half range resolution bin is considered. In GNSS-based passive radar, the best range resolution we can have is c/2B_w in the quasi-monostatic case, and the range-curvature term is $\frac{\lambda }{2}\left(f_{r0}+f_{r,v}\right)\cdot {\left(\frac{T_c}{2}\right)}^2$. So, this assumption determines the upper bound of coherent time T_c, which is:

$$T_c\le \sqrt{\frac{2c}{\lambda B_w\left(f_{r0}+f_{r,v}\right)}}$$

(29)

Based on the upper bound and lower bound of the coherent time, Figure 8a has shown the available coherent time for aircraft detection in GNSS-based passive radar, with the parameters listed in Table 1 and SNR_min equal to 12.8 dB. Then, with coherent time T_c equal to 5 s, the result of power-budget analysis is shown in Figure 8b. As shown, the detection range of aircraft using GNSS-based passive radar is larger than 80 km, and it will be increased twofold in the RCS enhancement area. Furthermore, with the advanced target detection method, like TBD, the detection range can easily reach 200 km, as shown in Figure 8b.

4.2.2. Searching Step of Doppler Rate

In (23), the residual quadratic phase error still exists because of the error of estimated Doppler rate $\Delta f_r$, and it mainly depends on the searching step of Doppler rate $\Delta {\widehat{f}}_r$. Generally, the smaller searching step we choose, the smaller the error of estimated Doppler rate. However, choosing a small searching step means the searching number is accordingly increasing, which will significantly increase the computational load. Therefore, choosing the parameter $\Delta {\widehat{f}}_r$ is very important in the proposed algorithm.

Due to the error $\Delta f_r$, the 2-D coherent integration gain is reduced, and the energy loss at the position of target motion parameters (f_d,v, R_ref) is:

$$E\left(\Delta f_r\right)=\frac{S_o\left(f_{d,v},R_{ref}\right)}{{S_o\left(f_{d,v},R_{ref}\right)|}_{\Delta f_r=0}}=\frac{\left|{\displaystyle \sum _{m=-N_A/2}^{N_A/2-1}\exp \left\{-j\pi \Delta f_r{\left(m{T}_p\right)}^2\right\}}\right|}{N_a}<\frac{\left|{\displaystyle \sum _{m=-N_A/2}^{N_A/2-1}\exp \left\{-j\pi \Delta {\widehat{f}}_r{\left(m{T}_p\right)}^2\right\}}\right|}{N_a}$$

(30)

Based on (30), with the parameters listed in Table 1, Figure 9 shows the impact of searching step on the energy loss $E\left(\Delta f_r\right)$ of the final processing results. In Figure 9, as the searching step increases, the energy loss will increase significantly. In our case (T_c is 5 s), the energy loss (less than 0.2 dB) could be ignored when $\Delta {\widehat{f}}_r$ is 0.04 Hz/s. Besides, the variation of f_r,v is smaller than 0.4 Hz/s, as shown in Figure 2b, and as a result, the maximum searching number is only 10, which is absolutely acceptable.

4.2.3. Doppler Ambiguity

According to the analysis in Section 2.1, the Doppler frequency f_d,v varies from −1500 Hz to 1500 Hz. Therefore, the range of searching Doppler frequency is [−1500, 1500] Hz, and the Doppler ambiguity is inevitable because f_p of the GNSS signal is 1000 Hz. Let the ambiguous Doppler frequency be:

$$f_{d,v}^a=f_{d,v}^{}+I\cdot f_p$$

(31)

where $f_{d,v}^a$ is the ambiguous Doppler frequency, and I is the ambiguous Doppler index. So, the range and Doppler profiles of the final coherent integration processing result can be written as:

$$\begin{array}{c}{S}_o\left(f_{d,v}^a,R\right)=A{Q}^C{\displaystyle \sum _{m=-N_A/2}^{N_A/2-1}D\left[R-R_{ref}-\delta \left(I\cdot f_p\right)+\lambda \cdot I\cdot m\right]}\\ S_o\left(f_v,R_{ref}\right)=A{Q}^C{\displaystyle \sum _{m=-N_A/2}^{N_A/2-1}D\left[-\delta \left(f_v-f_{d,v}\right)+\lambda \left(f_v-f_{d,v}\right)m{T}_p\right]}\cdot \exp \left\{j2\pi \left(f_v-f_{d,v}\right)m{T}_p\right\}\end{array}$$

(32)

where the residual quadratic phase term has been ignored. With the parameters listed in Table 1, the interpolated range and Doppler profiles are shown in Figure 10. As shown, energy loss exists in the range dimension, and the final result cannot be focused in the range-dimension, because of the ambiguous Doppler frequency. To evaluate the impact of Doppler ambiguity, the blind Doppler sidelobe ratio (BDSR) can be defined as:

$$BDSR\left(I\right)=\frac{S_o\left(f_{d,v}^a,R_{ref}\right)}{S_o\left(f_{d,v}^{},R_{ref}\right)}=\frac{{\displaystyle \sum _{m=-N_A/2}^{N_A/2-1}D\left[-\delta \left(I\cdot f_p\right)+\lambda \cdot I\cdot m\right]}}{N_A}$$

(33)

where I is equal to ±1, ±2. When the coherent time is 5 s, N_A = 5000, BDSR(±1) = −17.06 dB, and BDSR(±2) = −27.40 dB. As the BDSR is very low, the blind Doppler sidelobe will submerge in the noisy background. Therefore, both the range/Doppler profile and BDSR results demonstrate that the proposed algorithm has a good performance even in the Doppler ambiguity situation, and the problem of blind speed is also overcome.

4.2.4. Computation Analysis about the Azimuth Coherent Processing

In GNSS-based passive radar, the long coherent time is needed and the range Doppler searching is implemented. So, the azimuth coherent processing could be time-consuming. In this section, the computation complexity analysis is given.

In digital processing case, one complex multiplication needs six floating-point (FLOP) operations, and the signal FFT and IFFT operation with a length of N needs 5Nlog₂(N) FLOP. Therefore, based on (13) and (17), the azimuth processing computation load is:

$$C_1=N_K\times \left[6N_A{N}_R+36N_V{N}_A{N}_R^{}\right]$$

(34)

Then, the fast implementation method based on CZT for azimuth coherent processing is presented. According to (20), the fast azimuth processing computation load is:

$$C_2=N_K\times \left[6N_V{N}_R+6N_A{N}_R+15N_A{N}_R{\log }_2{N}_A\right]$$

(35)

Assume that the computational efficiency γ is defined as C₁/C₂. Normally, we let N_V = N_A, and with the condition of the same parameters as listed in Table 1, the computational efficiency γ is:

$$\gamma =\frac{C_1}{C_2}=\frac{6N_A{N}_R+36N_V{N}_A{N}_R^{}}{6N_V{N}_R+6N_A{N}_R+15N_A{N}_R{\log }_2{N}_A}\approx \frac{36N_A}{12+15{\log }_2{N}_A}\approx 900$$

(36)

Therefore, the computational complexity of the proposed method is significantly reduced by the CZT operation.

5. Conclusions

GNSS-based passive radar is a useful tool to detect aircrafts, because of its system advantages and signal characteristics. According to the power budget analysis, the detection range of aircrafts is ranging from 80 km to 150 km, and it can still be improved if advanced moving target detector is adopted. However, during the echo signal processing, there are great challenges to achieve coherent integration gain in both range and azimuth dimensions, which will limit the performance of GNSS-based passive radar significantly.

Therefore, a novel 2-D coherent integration processing and detecting algorithm is proposed for aircraft detection. In this paper, a modified RFT with range-walk removal and Doppler rate estimation is performed for azimuth long time coherent processing, and its fast implementation method based on the Chirp-Z transform is also presented. Furthermore, a modified matched filter with a shifting Doppler frequency in each Doppler line is performed for range compression. Hence, 2-D coherent integration gain is obtained, and the range-Doppler focused result has a sufficiently high SNR for detection. The excellent performance of the proposed algorithm is demonstrated by numerical simulation experiments. Moreover, as the experimental results shown, the blind Doppler sidelobes are weak enough to submerge in the noise background, which demonstrate that the proposed algorithm has good performance even when Doppler ambiguity exists. With the estimated motion parameters of the aircraft, the aircraft image could be obtained by GNSS-based SAR imaging formation algorithm. Although the fast implementation method based on CZT is used in the azimuth coherent processing, the computation load of the proposed algorithm is still large because of the long coherent time. Therefore, in real-data experiments, parallel computing based on FPGA or GPU is highly recommended.

Although GNSS signals give low power density at the aircraft targets, their abundance and the variety of geometric configurations from different satellites provide two major advantages: (1) The bistatic angle giving largest RCS can be selected to improve the radar detection range; (2) the multiple observation geometries allow the algorithm to obtain not only the range and Doppler parameters of aircraft but also its position and velocity. Therefore, based on experiments and power budget analysis, this paper has also highlighted the potential of GNSS-based passive radar as an assistive tool for air traffic control.

The near future work in this direction is to test the proposed algorithm with real GNSS-based passive radar echo signals. Later work will explore aircraft target detection and tracking in low SNR, as this could dramatically increase the GNSS-based passive radar’s range of coverage.