GNSS-Based Passive Inverse SAR Imaging

The utilization of global navigation satellite system (GNSS) signals for remote sensing has been a hot topic recently. In this article, the feasibility of the GNSS-based passive inverse synthetic aperture radar (P-ISAR) is analyzed. GNSS-based P-ISAR can generate the two-dimensional image of a moving target, providing an estimation of the target size, which is very important information in target recognition. An effective GNSS-based P-ISAR moving target imaging algorithm is proposed. First, a precise direct path interference (DPI) suppression method is derived to eliminate the DPI power in the detection channel. Then, the P-ISAR signal processing method is established. Due to the large synthetic aperture time, the Doppler profile of the ISAR image will defocus if directly performing the Fourier transform. As a solution, a parametric autofocusing and cross-range scaling algorithm is specially tailored for the GNSS-based P-ISAR. The proposed algorithm cannot only focus and scale the ISAR image, but also provide an estimation of the cross-range velocity of the target. Simulation with an airplane target is designed to test the signal processing method. Finally, an experiment is conducted with a civil airplane as the target and GPS satellites as the illumination source. Focused ISAR image is successfully acquired. The estimated length and velocity of the target are approximately consistent with ground truth, which are obtained by the flight record. The potential of the GNSS-based P-ISAR on multistatic operations is also illustrated by the fusion of the ISAR images obtained using different satellites as illumination sources.


I. INTRODUCTION
I N ADDITION to providing positioning and timing services, the global navigation satellite system (GNSS) also provides high-quality L-band radiation sources with real-time global coverage. GNSS signals are constantly beamed to Earth and reflected back into space with an imprint of the Earth's surface conditions. Recently, the secondary use of GNSS signals has been a hot topic in remote sensing community [1]. The main idea Manuscript  is to make use of GNSS reflected signals, which are considered as multipath interference in positioning application, to derive the geophysical parameters of the Earth's surface. Relevant topics include GNSS reflectometry (GNSS-R) [2], [3], [4], [5], [6], [7], GNSS-based synthetic aperture radar (SAR) [8], [9], [10], [11] and SAR Interferometry (InSAR) [12], [13], and GNSS-based passive radar [14], [15], [16], [17], [18], [19], [20], [21]. GNSS-R is the earliest and most mature technic related to the secondary use of GNSS signal. It utilizes the delay-Doppler map of GNSS reflected signals to retrieval the geophysical parameters of Earth's surface, such as ocean wind speed, ice-layer density, forest biomass, and ocean surface altitude. Several dedicated GNSS-R space missions have been successfully launched in the past few years [6], [7].
GNSS-based SAR and SAR Interferometry (InSAR) aims at obtaining the two-dimensional microwave image and digital elevation model of the interested area. In [15] and [16], the feasibility of using a GNSS-based passive radar to detect moving ship target is illustrated, and long-time integration algorithms for noncooperative target are proposed. Relevant works can also be found in [20] and [21], which show the possibility of detecting air and ground moving targets using a GNSS-based passive radar. In [18], a passive radar imaging algorithm for the ship target is proposed, and the potential of the GNSS-based PBR in moving target classification is illustrated.
Inverse SAR (ISAR) is a powerful radar signal processing that is able to form the imagery of noncooperative moving targets [22], [23]. Passive ISAR (P-ISAR) utilizes illuminators of opportunity as the radiation source, and has also becomes a research hotspot recently [24], [25], [26], [27]. In this article, the feasibility of the GNSS-based P-ISAR is analyzed and tested by experiments. Generally, radar detection can provide the decision of whether a target is present or not, and a rough estimation of its speed. SAR can provide the two-dimensional radar image of the target. However, for high-speed targets, its SAR image would defocus due to the mismatch of the imaging parameters. GNSS-based P-ISAR processing starts from the detection result of the GNSS-based passive radar, and goes further to obtain the focused image of the noncooperative moving target. It can extract additional information about the target, such as the cross-range length and a finer estimation of the velocity, which are very meaningful for target recognition.
Due to the removal of transmitting device, GNSS-based P-ISAR has several advantages shared with other passive radar, such as low-cost, low-power, light-weighted, no requirement for additional spectrum allocation, and covert operations. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ The property of the GNSS signals also enables the global and real-time working ability of the GNSS-based P-ISAR. Finally, GNSS timing service provides great convenience for bistatic synchronization.
The problems of the GNSS-based P-ISAR lie in two aspects. First, as continuous wave radar, GNSS-based P-ISAR suffers from severe direct path interference (DPI). The detection antenna unavoidably receives a strong directly received transmission. The target echo must compete with DPI, which experiences only transmission path attenuation. Due to the BPSK modulation of GNSS ranging code, its side-lobes do not attenuate. In some cases, even the side-lobe of the DPI can be stronger than the target echo. Therefore, a proper DPI suppression method must be derived for the GNSS-base P-ISAR. Second, due to the low effective isotropic radiated power (EIRP), the GNSS-based P-ISAR typically requires a long synthetic aperture time. As analyzed in the previous work [21], the required integration time for the GNSS-based passive radar can reach as long as 10 s depending on the detection range. Consequently, the Doppler profile of the ISAR image for high-speed targets will defocus. Generally, the range resolution of GNSS signals is fairly limited [28]. Hence, the precision of the Doppler phase obtained by envelope alignment is not sufficient for high-order phase compensation. Therefore, an effective autofocusing algorithm should be established for the GNSS-based P-ISAR.
In this article, the geometric and signal models of the GNSSbased P-ISAR are established. Based on the models, a signal processing algorithm is specially tailored for the GNSS-based P-ISAR. First, a precise DPI suppression method is applied to the detection channel signal, then the range compression and translational motion compensation (TMC) is performed. In order to obtain the final focused ISAR image, a parametric autofocusing and cross-range scaling algorithm is applied to the target echo. The cross-range length and velocity estimation of the target is obtained, which can help improve the target recognition capability of the GNSS-based passive radar. Simulation and experiment using an airplane as the target are conducted to confirm the feasibility of the proposed algorithm. The multistatic operation potential of the GNSS-based P-ISAR is also illustrated by the noncoherent fusion of the ISAR images obtained by using different satellites as the illumination source.
The rest of this article is organized as follows. The signal processing method is proposed in Section II. The simulation and experiment are presented in Section III. Finally, Section IV concludes this article.

A. Power Analysis of DPI
In this section, the relative power level of the DPI and target return is analyzed. According to the bistatic radar equation [29], the power of the DPI received by the detection channel antenna can be expressed as  where P t , G t , and R t are the transmit power, the transmit antenna gain, and the transmit range, respectively, λ is the carrier wavelength, EIRP is the EIRP of GNSS. G rd is the receiving antenna gain at the angle of arrival (AOA) of the DPI.
The power of the target echo received by the detection channel antenna can be expressed as where σ is the radar cross section (RCS) of the target, G r is the receiving antenna gain at the AOA of the target, R r is the receiving range.
Comparing (1) and (2), the ratio between the power of the DPI and the target echo in the detection channel is Generally, the main-lobe of the detection antenna is pointed toward the target, and the AOA of the DPI usually falls into the side-lobe, which means G rd ≤ G r . Fig. 1 shows the curve of the power ratio between the DPI and the target return as a function of detection range, where the utilized parameters are listed in Table I. As can be seen, the DPI collected by the detection channel is much stronger than the target echo. Fig. 2 shows the ambiguity function of the GPS L1 ranging code. The delay and Doppler profiles are plotted in Fig. 3. As can be seen, due to the BPSK modulation, the side-lobes of the Doppler cut do not attenuate. It can be calculated that the peak-to-side-lobe ratio of the autocorrelation function for GPS L1 ranging code is 24 dB. When the power of the DPI is more than 24 dB higher than that of the target echo (indicated by the red dashed line in Fig. 1), even the side-lobes of the DPI can be stronger than the target echo. Therefore, an effective  DPI suppression method must be developed for the GNSS-based P-ISAR.

B. Signal Processing Method
The system geometry of the GNSS-based P-ISAR is shown in Fig. 4. The receiver is fixed on the ground. Two receiving channels are exploited, i.e., the reference channel for bistatic synchronization, and the detection channel for target imaging. The GNSS satellite continuously transmits BPSK modulated signals, which are reflected by the moving target and collected by the antenna of the detection channel. Different from GNSSbased passive SAR, which utilizes the motion of the receiver to form a synthetic aperture, GNSS-based P-ISAR relies on the motion of the target to realize cross-range resolution.
The transmitted signal of the GNSS satellite can be expressed as where t represents time, and f c represents the carrier frequency. D(·) represents GNSS data code. C(·) represents GNSS ranging code. N n=0 rect( represents the ranging code pulse train with period T p and pulse number N. The amplitude and initial carrier phase has been omitted since they do not have impact on the analysis.
The signal received by the reference channel is the delay of the GNSS transmitted signal, which after carrier demodulation can be expressed as where R d (t) is the propagation range from the GNSS satellite to the reference antenna. The point target echo received by the detection channel can be expressed as where R(t) represents the total propagation range of the target echo. R t (t) is the range from the transmitter to the target, and R r (t) is the range from the target to the detection antenna. The block diagram of the proposed GNSS-based P-ISAR imaging algorithm is shown in Fig. 5. As analyzed in Section II-A, the DPI and even its side-lobes are much stronger than the target echo. Consequently, a DPI suppression method is first performed to the detection channel signal.
1) DPI Suppression: DPI suppression is an important research topic of the bistatic continuous wave radar. The most utilized DPI suppression methods include spatial filtering method, the CLEAN type algorithms [30], [31], and the extensive cancellation algorithm (ECA) [32]. The spatial filtering utilizing the adaptive beam forming technique to steer zeros of the antenna pattern in specific directions to cancel interference sources. However, the implementation requires a phased array antenna, and the suppression performance may be not sufficient for strong interference. The CLEAN method is first proposed in radio-astronomy to detect the weak signal superimposed by the strong interference [30]. Recently, it is introduced in radar signal processing to improve the performance in image dynamic range [31]. However, the derivation of the weighted correlation coefficient for the CLEAN method requires the analytical expression of the reference signal, which is hard to obtain for the GNSS signals. The ECA projects the reflect channel signal into a subspace orthogonal to the extensive interferences, and are widely used for the good cancellation performance and convergence. However, the performance of the ECA decrease for timevarying signals. In this article, a DPI suppression method based on the high-precision interference reconstruction is proposed for the GNSS-based P-ISAR, which can also perform well when the signal amplitude fluctuates with time.
The proposed DPI suppression algorithm is based on the fact that the direct path signals in the reference channel and the detection channel are highly correlated, as long as the detection channel antenna is located close to the reference channel antenna. This condition is very similar to the InSAR system [33], which utilizes the interference property of the SAR images generated by two closely located antennas. The DPI in the detection channel can be expressed as where k and Δφ are the relative amplitude and phase errors (APE) between the direct path signals in two channels. If the direct path signal in the reference channel can be replicated, and the APE can be accurately estimated, the DPI in the detection channel can also be reproduced and removed precisely. The block diagram of the proposed DPI suppression method is shown in Fig. 6. It can be divided into two steps, i.e., the direct path signal replication step, and the APE estimation step. Generally, GNSS signal has a three-layer modulation structure [28]: the bottom layer carrier phase modulation, the middle layer ranging code modulation, and the top layer data code modulation. In order to replicate the direct path signal, all of the three modulations should be precisely reconstructed.
The direct path signal in the reference channel is replicated aided by GNSS tracking measurements, as shown in Fig. 6. First, the carrier phase modulation is generated by the integration of the Doppler frequency measurement. In order to eliminate the impact of system noise, the tracked carrier frequency is first low-pass filtered, which is accomplished by a polynomial fitting. Then, the pseudorange measurement is utilized to calculate the instantaneous code phase. Similarly, the tracked pseudorange is low-passed filtered by polynomial fitting to smooth the noise. Once the instantaneous code phase is obtained, the ranging code modulation can be reconstructed. Finally, the data code modulation is generated according to the prompt code correlation result. The prompt code correlation result contains not only the data code modulation message, but also the instantaneous amplitude of the signal, and consequently, the amplitude of the direct path  signal can also be retrieved in this step. The denosing of the prompt code correlation measurement is accomplished by an innerchip mean filtering. Fig. 7 shows an example of the Doppler frequency measurement with a length of 100 s, which is plotted with a black solid line. The filtered result is also plotted in the figure with a red dashed line. Since the Doppler frequency of the reference channel signal is totally caused by satellite motion, a polynomial fitting can well recover the true value from the Noise-contaminated measurement. Fig. 8 shows an example of the prompt code correlation measurement, which is plotted with a black solid line. The sign of the correlation value indicates the phase of the data code, whereas the jitter of the amplitude represents the noise contamination. The filtered result is also plotted in Fig. 8 with a red dashed line. Between two sign transitions, the correlation value is mean filtered to smooth the noise.
DUE to the different transmit path and channel characteristics, the generated direct path signal has an APE relative to the DPI. The precise estimation of the APE is the key step for effective DPI suppression. The block diagram of the proposed APE estimation method is shown in Fig. 9. The main idea is to coherently integrate a certain number of pulses for the two channels, respectively, and then calculate the ratio between the  integration results. The complex valued ratio is chosen as the estimation of the APE. The multipulse integration can reduce the impact of the noise on the APE estimation, and the number of the integrated pulses can be set according to the desired estimation precision and processing efficiency. The direct path signal replica obtained in previous step is multiplied by the APE estimation to generate the final DPI estimation. Fig. 10 shows an example of the final estimated DPI. The DPI suppression can be accomplished by subtracting the estimated DPI from the detection channel signal.
2) ISAR Processing: After DPI suppression, the majority of the DPI power in detection channel is eliminated. The P-ISAR imaging can then be performed. The data code is first demodulated, and then the one-dimensional data is rearranged to form the two-dimensional radar matrix. According to [21], the radar matrix at this step can be expressed as where η represents the slow-time, τ = t − nT p is the fast-time, T a = NT p is the total length of the signal being processed.
In (8), the Doppler phase modulation within the ranging code is omitted, i.e., the stop-and-hop model is assumed. However, according to Fig. 2, the autocorrelation peak of the GNSS ranging code drops rapidly with Doppler frequency mismatch. Generally, the Doppler frequency caused by satellite motion can reach as large as ±6 kHz, much higher than the first null frequency of the ambiguity diagram. Thus, the satellite motion must be compensated during range compression. The range compression filter is designed as where f dc is the central Doppler frequency of the satellite relative to the receiver. The range compression is accomplished by the correlation between the range pulses and the range compression filter where represents the correlation operator. P (·) represents the autocorrelation function of the GNSS ranging code. At this step, the satellite motion has been compensated, and the synthetic aperture is mainly caused by the motion of the target. Generally, the motion of the target can be divided into two parts, i.e., the translational motion and the rotational motion [22]. The rotational motion generates the synthetic aperture angle, which is the source of the azimuth resolution capability of ISAR system. The translational motion causes range walk and phase errors of the target echo, which should be properly compensated. The TMC consist of a range alignment and a phase adjustment [22]. After TMC operation, the target motion can be equalized to the turntable model, and the signal can be expressed as where R rot (η) = R res (η) − R TM (η), and R TM (η) is the range caused by the translational motion of the target. Under the turntable model, the target can be viewed as rotating around the central axis. The Doppler frequency of a scatter is proportional to the cross-range of that scatter relative to the central axis. The two-dimensional ISAR image of the target can be obtained through a Fourier transform along the crossrange dimension. However, in the GNSS-based P-ISAR, the synthetic aperture time is typically very long to enhance the signal-to-noise ratio (SNR) of the ISAR image. The acceleration of the rotational angular velocity of the target will cause Doppler spectrum defocus, and hence the degradation of the ISAR image. Even if the rotational angular velocity maintains constant, the Doppler spectrum will still defocus with a long synthetic aperture time. In order to improve the image quality, the nonlinear Doppler phase caused by target motion must be properly compensated. Since the resolution of the GNSS ranging code is fairly poor (tens of meters to hundreds of meters according to the bandwidth of the utilized signal), the Doppler phase retrieved by envelop tracking is not precise enough to achieve this purpose. The precision requirement on range measurement for phase compensation purpose is related to the signal wavelength, much stricter than that for range alignment purpose.
In order to solve this problem, a parametric autofocusing and cross-range scaling algorithm is proposed in this article. The parametric autofocusing step generates the focused Range-Doppler image of the target, and obtains an estimation of the cross-range velocity of the target. The cross-range scaling step scales the Doppler frequency to the cross-range utilizing the estimated velocity acquired in the previous step, and generates the final ISAR image.
Expanding R rot (η) at η = 0 and keeping only those up to the second order, we have where R 0 = R rot (η)| η=0 , f d and f r are the equivalent Doppler frequency and Doppler rate caused by the rotational motion, respectively. The signal after TMC in (13) can then be expressed as The linear and second-order term in signal envelope is neglected since they are typically far smaller than the range resolution. The Doppler frequency f d indicates the target cross-range position after Fourier transform, whereas the Doppler rate f r causes Doppler spectrum defocusing, which should be compensated. Since the satellite motion has been compensated yet, the Doppler rate f r is only related to the motion of the target. For the GNSS-based P-ISAR, since the satellite orbit is very high (MEO and GEO orbit), the Doppler rate f r is mainly caused by the target motion relative to the receiver. It can be calculated that where v ⊥ is the tangential velocity of the target relative to the receiver (as marked in Fig. 4), R r0 is the range from the target to the receiver at the central synthetic aperture time. A bistatic range scaling procedure is needed to scale the bistatic range delay R 0 to the ground range R r0 . Considering the specific detection configuration of the GNSS-based P-ISAR, two solutions are provided here. The first is a single-satellite solution which requires a narrow beam receiving antenna. The single-satellite solution is based on the fact that the transmit range is far larger than the receiving range for the GNSS-based P-ISAR. Thus, the line of sight from the transmitter to the target is nearly parallel to that from the transmitter to the receiver, as shown in Fig. 11. Letting point P be foot of perpendicular from the receiver to the transmitting line of sight, we have where R s1 is the range from the transmitter to point P, R s2 is the range from point P to the target, β is the bistatic angle, which can be calculated by where α and θ are the elevation angle and horizontal angle of the GNSS satellite, respectively, which can be determined by the GNSS positioning service, ϕ is the AOA of the target, which can be determined by the pointing direction of the receiving antenna. The bistatic scaling can be accomplished utilizing (18) to (20). The precision of the single-satellite solution is limited to the angular resolution of the receiving antenna, which is typically fairly rough. The second approach is a multistatic solution. Generally, if more than 3 GNSS satellites can be utilized to detect the target, the three-dimensional position of the target can be estimated by solving the range delay equation set. This multilateration technique can be treated as a modification of the standard GNSS positioning procedure, which has already been discussed and illustrated by experiments in [21] for the GNSS-based passive radar. Once the target position is obtained, the bistatic angle β can be calculated, and the bistatic range scaling can be performed utilizing (18). The precision of the multistatic solution is limited only to the ranging precision and the spatial distribution of the utilized satellites. Since GNSS is dedicatedly designed for positioning, the multistatic solution is typically more precise than the single-satellite solution.
Since the accurate velocity of the noncooperative target is unknown, a parametric autofocusing method is introduced to eliminate the Doppler spectrum defocusing. The parameter needed to be tuned is the tangential velocity of the target relative to the receiver, and the phase compensation function is The discrete form Doppler profile after phase compensation is (22) The Doppler profile contrast function is chosen as the objective function that needs to be maximized, which is defined as . (23) The search region of the tangential velocity where v max is the maximum velocity estimation of the detected target. The search step is limited by the synthetic aperture time The search parameter v * ⊥ corresponding to the maximized objective function is chosen as the estimation of the tangential velocity of the target. The resulted Range-Doppler image of the target with maximum contrast is where AFFT[·] represents the fast Fourier transform in the crossrange direction. Once the tangential velocity and Doppler rate of the target has been estimated, the cross-range scaling of the ISAR image can be accomplished. For a scatter with a Doppler frequency f η , its corresponding cross-range position can be determined as . (27) The estimation of the rotation velocity (or the tangential velocity v ⊥ ) is the key for the cross-range scaling of the ISAR image. The conventional way to estimate the target velocity is by the envelope tracking of the range compressed echoes [34], [35]. Once the range migration of the target is obtained, the velocity of the target can be retrieved. This method works well for the wideband ISAR system thanks to the high range resolution. However, the range resolution of the GNSS ranging code is fairly poor due to the limited bandwidth, the velocity retrieved by envelop tracking is typically not sufficient for precise scaling of the ISAR image. For the proposed algorithm, the tangential velocity of the target is estimated by the phase information of the echoes, which is typically more precise than the envelopbased estimation method. Supposing the phase error caused by Doppler rate estimation error is Δφ, The corresponding Doppler rate error is where T a is the integration length. The relation between the Doppler rate and the tangential velocity is The velocity estimation error can be calculated by For a target with a cross-range position x, the scaling error caused by the velocity estimation error can be determined by

III. SIMULATION AND EXPERIMENT
In this section, first, a simulation is performed with an airplane model as the input to demonstrate the signal processing method. Second, an experiment is conducted with a civil airplane as the target to further confirm the feasibility of the proposed GNSSbased P-ISAR algorithm.

A. Simulation
The simulation input is a civil airplane model, as is shown in Fig. 11. The length of the airplane is set to be 46 m, and the width is set to be 40 m. The positions of the GNSS satellite, the target and the receiver are listed in Table II. Their velocities are listed in Table III. Some parameters used in the experiment are listed in Table IV.   TABLE II  POSITIONS OF THE TRANSMITTER, THE RECEIVER, AND THE TARGET   TABLE III  VELOCITIES OF THE TRANSMITTER, THE RECEIVER, AND THE TARGET   TABLE IV  The range compressed result is shown in Fig. 12. Since the satellite motion has been compensated, the range migration is mainly caused by motion of the target. According to Tables II and III, the central Doppler frequency of the target is expected to be fairly small under the simulation parameter settings. Due to the interference of the targets' echo, amplitude fluctuation can be observed for the range compressed signal appear in the same range cell. TMC is then performed, which consist of a range alignment and phase adjustment. Fig. 13 shows the range Doppler image of directly performing a Fourier transform in the cross-range direction (the amplitude has been normalized to the mean noise power level, and is shown in dB scale). The Doppler dimension has been scaled to crossrange by (27) using the estimated v * ⊥ and f * r in the following step.  Due to the long integration time and large tangential velocity of the target relative to the receiver, significant defocusing can be observed in the Doppler profile of Fig. 13. the Doppler profile defocusing not only will decrease the SNR of the target, but also will affect the determination of the target size. The crossrange profile of the central range cell is plotted in Fig. 14 (black solid line). As can be seen, the length of the target image has expanded to more than 300 m, which will seriously interfere with the identification of the target. The proposed parametric autofocusing algorithm is then performed to the echo signal. The Doppler profile contrast of the central range cell as a function of the tuning parameter, the tangential velocity of the target relative to the receiver v ⊥ , is plotted in Fig. 15. As can be seen, the Doppler profile contrast function has a peak value when v ⊥ = 195 m/s, and this value is chosen as the optimum estimation of the tangential velocity of the target, i.e., v * ⊥ = 195 m/s. The single-satellite solution is utilized in the ground range estimation. Since no pointing error is added in the simulation, the estimated ground range is well consistent with the true value. The optimum estimation of the Doppler can also be acquired by (16), which is f * r = 29.8 Hz/s. The estimated parameters are summarized in Table V. The ISAR scaling can then be accomplished by (27) using these estimated parameters.
The final scaled ISAR image of the target is shown in Fig. 16. As a comparison, the Doppler profile of the autofocused ISAR image is also plotted in Fig. 14. As can be seen, the Doppler profile is well focused. The estimated target size after autofocusing as well as their relative error to the ground truth are listed in Table VI. As can be seen, both the estimated range and cross-range lengths are approximately consistent with the true values. The relative error may be caused by the spectral leakage due to signal truncations. In the simulation, the target is set to have a tangential velocity, which contributes only to the cross-range resolution. For the range dimension, the resolution is achieved by measuring the transmission delay of the reflected pulses. The focused ISAR image is first up-sampled in both dimension, and then the target size is estimated by counting the pixels occupied by the image.

B. Experiments
Experiments are conducted to further confirm the feasibility of the proposed GNSS-based P-ISAR imaging algorithm. The experiment is conducted near the Nan-Yuan airport, Beijing, China, in January 5, 2019. An illustration of the experimental geometry is shown in Fig. 17. The target is a landing airplane which is fly northward. The receiver is located at the east of the landing path. The line of sight of the detection channel antenna is nearly perpendicular to the landing path. The distance from the receiver to the target is about 300 m. Under the experimental scenario, the long side of the airplane is approximately parallel to the cross-range direction of the ISAR image. Consequently, the cross-range length of the ISAR image is an estimation of the airplane length. Fig. 18 shows photo of the experimental setups and the target token during the experiment. The reference channel antenna is positioned on the ground, and the detection channel antenna is mounted on a tripod a few meters away from the reference antenna.
The data are collected when the target passes the main-lobe of the detection antenna. The airplane is successfully detected whose flight information can be found in the flightradar24 website. According to the flight information, the aircraft model is Airbus A321 whose length is 44.5 m. The ground velocity of the airplane at the detection moment is 66 m/s. Fig. 19 shows the sky plot of the visible GNSS satellites acquired by  the reference channel signal. The sky plot is plotted in the local geodetic east-north-up coordinate system. The spokes indicate the horizontal angle of the satellite, whereas the circles indicate the elevation angle.
Limited by the physical size, the exploited detection antennas cannot determine the precise pointing direction (more than 20°b eam width in horizontal direction). Therefore, the multistatic solution for bistatic range scaling is performed. In the reported experiment, two detection antennas are exploited (as is shown in Fig. 19   target is estimated by solving the ranging equation set for the 4 GNSS satellites utilizing the method proposed in [21].
The satellites utilized for ISAR imaging are marked with red circles in Fig. 20, which are GPS SVN 25 and SVN 32. Both of these two satellites belongs to the generation Block IIF satellites, transmitting the newly added GPS L5 signal that has a better range resolution than the L1 signal. Fig. 21 shows the range compressed result of the detection channel signal utilizing SVN 25 satellite as the illumination source. The length of the presented signal is 3 s. As can be seen, the DPI appearing at zero delay region is stronger than the target return. The power of the DPI is relatively stable within the total length of 3 s. The target echo appears at the moment of about 1 s and disappears at about 2 s. The proposed DPI suppression algorithm is then performed. Fig. 22 shows the range compressed result after DPI suppression. As can be seen, the signal power  of the DPI has been greatly suppressed, with the power of the target echo being nearly not affected. Since the satellite motion has been compensated in advance, and the velocity of the target is nearly perpendicular to the radar line-of-sight, the target echo in Figs. 21 and 22 shows no significant range migration in the detection duration. Due to the short detection range and the large tangential velocity of the target, the synthetic aperture time is set to be less than 3s. When the detection range increases, the required synthetic aperture time will also increase.
The TMC is performed to compensate the translational motion of the target, which includes range alignment and phase adjustment. After TMC, target echo envelopes are precisely aligned, and the Doppler phase of the echo caused by translational motion is compensated. When the synthetic aperture time is short, the ISAR image of the target can be obtained by directly performing an azimuth Fourier transform to the target echo after TMC. However, for the GNSS-based P-ISAR, in order to enhance the target SNR, the synthetic aperture time is relatively long. The Doppler spectrum of the target echo will defocus when directly performing the azimuth Fourier transform. Fig. 23 shows the  ISAR image of the airplane target, which is obtained by directly performing an azimuth Fourier transform, and scaled utilizing the parameters estimated in the following steps (the amplitude has been normalized to the mean noise power level, and is shown in dB scale). As can be seen, the defocusing of the spectrum in cross-range direction is very severe, and the estimated length of the airplane has been extended to about 80 m, which is nearly double the true length of 44.5 m.
The parametric autofocusing method is then performed to the signal after TMC. Fig. 24 shows the Doppler profile contrast as a function of the tuning parameter, the tangential velocity v ⊥ of the target relative to the receiver. As can be seen, the Doppler profile generates a peak value at v ⊥ = 59.3 m/s, which is chosen as the optimum velocity estimation v * ⊥ . The estimated range of the target relative to the receiver is R r0 = 280 m. The Doppler phase compensation function is established according to (21) using the estimated v * ⊥ and R r0 . The ISAR image of the airplane after parametric autofocusing is shown in Fig. 25. The Doppler frequency has been scaled to the cross-range according to (27). It can be seen that the Doppler spectrum of the target has been well  focused after autofocusing. The comparison between the crossrange profile of the ISAR image before and after autofocusing is plotted in Fig. 26. The estimated cross-range length of the target is 38.5 m, as listed in Table VII, approximately consistent with the true length of the airplane. As listed in Table VII, the estimation length is 6 m shorter than the true value. This may be because the flight path of the aircraft and the radar line of sight are not strictly vertical at the detection moment. Besides, the ISAR image indicates the backscattering intensity of the target, and for some structures of the airplane, such as the nose, the backscattering may be not strong under the experimental configuration. This can also result in the estimation error of the target length.
One of the advantages of the GNSS-based passive radar is the potential of multistatic operations. Fig. 27 shows the ISAR imaging results of the target airplane using GPS satellite SVN32 as the illumination source. As can be seen, the position of the strong scattering points is a little different from the result of using SNV25 signal. Different illumination angles may help to obtain more comprehensive information on the target. As an illustration, Fig. 28 shows the noncoherent image fusion result of the two ISAR images. As can be observed, the SNR of the image has a slight improvement. The estimated cross-range becomes 40.9 m, also more precise than the estimation obtained by the single source image.
High-precision registration is the premise of image fusion. The registration of the two ISAR images in range dimension is accomplished by the bistatic range scaling. The registration  in Doppler dimension is accomplished by a Doppler shift of the target echo according to their detected Doppler centroid frequency. The registered images are then normalized to their mean noise power level, and added together to form the fused image.

IV. CONCLUSION
In this article, the feasibility of the GNSS-based P-ISAR for moving target imaging is analyzed. An effective signal processing algorithm is proposed, which starts with a DPI suppression step to eliminate the DPI power in the detection channel signal. In order to eliminate the Doppler profile defocusing caused by the long synthetic aperture time, a parametric autofocusing and cross-range scaling algorithm is proposed. The proposed algorithm cannot only focus and scale the ISAR image, but also provide an estimation of the cross-range velocity of the target. The proposed GNSS-based P-ISAR algorithm can obtains the cross-range length and velocity estimations of the target, which are very valuable information for target classification, and will help improve the target recognition capability of the GNSSbased passive radar. The proposed DPI suppression and autofocusing algorithm can be conveniently adapted to the passive radar of using other illuminators of opportunity.