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

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.


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]. Based on the system geometry, the range history is the sum of the instantaneous receiver slant range RR(η) and transmitter slant range RT(η), which can be rewritten by Taylor series expansion.
( ) ( ) ( ) 2 2 λ η η η λ η η where η, Rref, λ, fd, and fr denote the illumination time, reference range, signal wavelength, target Doppler centroid, and Doppler modulated rate, and fd and fr are derived as follows [36]: , (2) where η = 0 means the illumination centroid time, fd0 and fr0 represent the Doppler centroid frequency and Doppler frequency modulated rate of the stationary target, and fd,v and fr,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 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.
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 r0 + f r,v (2) 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 fd,v and fr,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 fd,v and fr,v, using the parameters listed in Table 1. As shown in Figure 2, both fd,v and fr,v are greatly influenced by passive radar geometry and target motion direction. fd,v is mainly determined by the equivalent range velocity, and fr,v is mainly influenced by the equivalent azimuth velocity [36]. Normally, as depicted in Figure 2a, the range of fd,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), fp, 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, fr,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.

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: where superscript i represents using the i-th visible satellite as the transmitter, σ is the target radar cross-section, AT is the signal amplitude, Q(·) is the modulation code, and Tp, c, and λ denote signal period, speed of light, and signal wavelength, respectively. Then, set η = nTp + τ, tn = nTp, and the range history can be written as: 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).

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: 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: 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.
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: 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.

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.

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 0 t + r 0 , where r 0 is the initial target position and v 0 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]: As Doppler frequency is f v = 2v/λ, the range history can be rewritten by R(t) = λf v0 t/2 + r 0 , where f v0 is the target Doppler frequency. So, the RFT in the range-Doppler plane is: 2r Then, by defining the discrete searching range and Doppler parameters, the discrete form of traditional RFT has also been defined as [18]: 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 0 , v 0 ) or (r 0 , 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.

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. 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, fv). When the parameterized line (r0, v0) or (r0, fv0), 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.

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.

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:

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: where Q f (·) is the FFT form of the GNSS modulation code,A = σA T exp j 2πR re f λ . Then, the phase term and range-walk caused by f d0 and f r0 are compensated with the 2-D filter H p (t, f τ ).
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:

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: where • The first phase term in the first line of Equation (13) is QPE compensation, andf 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:f wheref 0 r , ∆f r and N k denote the initial Doppler rate, Doppler rate searching step, and Doppler rate searching number, respectively. Finally, whenf k r matches the true value of target Doppler rate or ∆ 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. 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 fd,v, the estimated Doppler rate equals fr,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(NANRNV), where NA, NR, and NV 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: (15) where Δ v f and τ Δf 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: where I and Nv denote the ambiguous Doppler index and non-ambiguous Doppler searching number. Based on this, Equation (13) can be rewritten as: 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 A N R N 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: where ∆ f v and ∆ f τ 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: 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: (17) Remote Sens. 2018, 10, 1164 Then, using the principle of the CZT, we have: where ⊗ denotes the convolution operation. Besides, based on the relationship between convolution and FFT, Equation (19) can be performed in digital processing as: 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 R N 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.
where ⊗ denotes the convolution operation. Besides, based on the relationship between convolution and FFT, Equation (19) can be performed in digital processing as: where subscript m means the operation in the azimuth direction, DFTm(·) and IDFTm(·) 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(NRNV), 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.

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 fd,v. This will introduce a time-varying phase in the range direction and change the response of range compression. The traditional cross-ambiguity function

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.
where subscript τ means the FFT operation in the fast time direction. By multiplying the modified matched filter H r ( f v , f τ ) and performing range IFFT operation, the signal after 2-D coherent processing is: where the superscript * represents the complex conjugate of H r ( f v , f τ ). Then, let R = cτ, and replacing the azimuth time integral with summation, the final result can be expressed as: 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.
is the factor caused by mismatched Doppler frequency in range compression, and δ(0) = 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 ∆ f r is zero, the signal S o f d,v , R re f is: 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 f d,v , R re f , f r,v . Therefore, by comparing the peak value of the processing result with different estimated Doppler rate,f r , the coherent integration gain can be achieved after the 2-D coherent integration processing.

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 −6 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.

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.

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. 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, 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. and the accuracy is better than one Doppler searching or range bin. Thus, the proposed method has a good performance on moving target focusing.

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.

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: , , , , , The connection between a target's location, velocity, and range Doppler parameters is: 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, all targets' locations and velocities are estimated accurately.

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: 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 n L 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: where T c is the coherent time, T p B w and f p T 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: 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 . So, this assumption determines the upper bound of coherent time T c , which is: 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. in Table 1 and SNRmin equal to 12.8 dB. Then, with coherent time Tc equal to 5 s, the result of powerbudget 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.
Based on (30), with the parameters listed in Table 1, Figure 9 shows the impact of searching step on the energy loss ( ) Δ r E f of the final processing results. In Figure 9, as the searching step increases, the energy loss will increase significantly. In our case (Tc is 5 s), the energy loss (less than 0.2 dB) could be ignored when Δˆr f is 0.04 Hz/s. Besides, the variation of fr,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.

Searching Step of Doppler Rate
In (23), the residual quadratic phase error still exists because of the error of estimated Doppler rate ∆ f r , and it mainly depends on the searching step of Doppler rate ∆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 ∆f r is very important in the proposed algorithm.
Due to the error ∆ 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: Based on (30), with the parameters listed in Table 1, Figure 9 shows the impact of searching step on the energy loss E(∆ f r ) 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 ∆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. Remote Sens. 2018, 10, x FOR PEER REVIEW 16 of 21

Doppler Ambiguity
According to the analysis in Section 2.1, the Doppler frequency fd,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 fp of the GNSS signal is 1000 Hz. Let the ambiguous Doppler frequency be: where , 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: where I is equal to ±1, ±2. When the coherent time is 5 s, NA = 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.

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: where f a d,v 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: 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: 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.

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 5Nlog2(N) FLOP. Therefore, based on (13) and (17), the azimuth processing computation load is: Then, the fast implementation method based on CZT for azimuth coherent processing is presented. According to (20), the fast azimuth processing computation load is: Assume that the computational efficiency γ is defined as C1/C2. Normally, we let NV = NA, and with the condition of the same parameters as listed in Table 1, the computational efficiency γ is: Therefore, the computational complexity of the proposed method is significantly reduced by the CZT operation.

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 GNSSbased 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

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 2 (N) FLOP. Therefore, based on (13) and (17), the azimuth processing computation load is: Then, the fast implementation method based on CZT for azimuth coherent processing is presented. According to (20), the fast azimuth processing computation load is: Assume that the computational efficiency γ is defined as C 1 /C 2 . Normally, we let N V = N A , and with the condition of the same parameters as listed in Table 1, the computational efficiency γ is: Therefore, the computational complexity of the proposed method is significantly reduced by the CZT operation.

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.
Doppler frequency is from −1500 Hz to 1500 Hz (shown in Figure 2), and the big Doppler shift is obviously not negligible, which means the Doppler frequency should be considered during the range compression processing.

Appendix B
The derivation of (13) is as follows: