An Image Formation Algorithm for Bistatic SAR Using GNSS Signal With Improved Range Resolution

Passive synthetic aperture radar (SAR) using Global Navigation Satellite System (GNSS) signals as the illuminators of opportunity has a relative coarse range resolution due to the narrow bandwidth of the navigation signals and the bi-static acquisition geometry, limiting its application for high-resolution requirement. In this paper, an image formation algorithm for range resolution improvement is investigated. A modified second order differentiation based method is applied on the range compressed signal of reference channel and surveillance channel, improving range resolution and keeping side-lobes at a lower level simultaneously, and the phase of processed signal is preserved by phase multiplication compensation for the following azimuth compression processing. To validate the effectiveness and benefits of the proposed algorithm in terms of range resolution and targets side-lobes, point targets simulation under different situations and real field experiment with strong scatters are conducted and discussed. Compared with the traditional GNSS-based passive SAR imaging algorithm, the results of the proposed method show that the imaging performance can be improved for better range resolution and lower side-lobes, benefitting for better distinguishing nearby targets and recognizing weak targets.


I. INTRODUCTION
Over the last few years, passive bistatic SAR (Bi-SAR) has attracted much attention from the research community. The separation of passive Bi-SAR system transmitter and receiver yields several configuration using different types of transmitters or receivers, including spaceborne, airborne, ground-based system [1]- [3]. Many different types of available illuminators of opportunity can be utilized as the transmission, such as digital audio broadcasting(DAB) [4], Digital Video Broadcasting Terrestrial(DVB-T) [5]- [10], WiFi [11], DTV [12], [13], and Global Navigation Satellite Systems (GNSS) [14]- [21]. The advantages of passive bistatic radar include low cost, free license, silence working and environmental friendly, making it feasible in the earth remote sensing fields.
The associate editor coordinating the review of this manuscript and approving it for publication was Junjie Wu.
Among all the illuminators of opportunity for passive bistatic radar, because of the abundant resources and global illumination of GNSS constellations, GNSS-based passive radar possesses a better ground coverage, and achieves global continuous observation, including the Polar Regions. The area can be illuminated from different perspective at the same time (typically six to eight satellites in a signal constellation simultaneously), and radar information space can be enhanced due to the different frequency band of satellites signals. Benefitting for these advantages, GNSS-based radar has been investigated and applied in many research fields, including surface change monitoring [22], ocean height measurement [23], oil slick detection [24], soil moisture retrieval [25], biomass estimation [26] and moving target detection [27]- [29].
The drawbacks of GNSS-based passive radar are low power budget and coarse resolution capability, because all GNSS satellites were not initially configured from earth VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ remote sensing. The weak signal from receiver can be enhanced by using long dwell time on the area. Long time coherent accumulation guarantees higher signal power in the processing, and improves the azimuth resolution by synthetic aperture radar technology simultaneously. However, the range resolution of GNSS-based passive SAR is limited by the transmitted signal bandwidth and the bistatic acquisition geometry. Although the best range resolution can be achieved in a quasi-monostatic configuration, for instance, the range resolution of conventional GPS-based SAR is about 150m for C/A code, which is not acceptable in most remote sensing scenarios. The research about range resolution improvement of passive radar mainly includes two strategies. The one is from the point of multistatic configuration systems, improving the azimuth and range resolution simulataneously. In 2015, Santi et al. [30] analyzed the point spread function of GNSS-based multistatic SAR theoretically and experimentally. Multi-pairs of bistatic images are incoherently combined, and the resolution can improved. In the same year, Zeng et al. [31] investigated the multi-angle imaging of BeiDou-2 based bistatic SAR. A region-based fusion method was applied on the single bistatic images, and higher quality fusion image was achieved, enhancing the image information space. In 2016, Santi et al. [32] introduced two types of fusion methods for the spatial resolution improvement of multistatic passive SAR. Extracting targets information using CLEAN technique before fusion processing and fusing multistatic images before targets information extraction were both testified, resulting in enhanced image information space and improved image performance respectively. In 2018, Nithirochananont et al. [33] investigated multistatic GNSS-based imaging using two different constellations of GNSS satellites, GPS and Galileo. Non-coherent addition was applied on the bistatic images, and targets geometric features were revealed including edge, shape and dimensions. Measwhile, object scattering variations can be exploited via a series of bistatic images. In 2020, Zheng et al. [34] investigated a multiple GNSS satellites coherently integrating imaging scheme, and range coordinate alignments were applied to perform azimuth coherent accumulation of multiple satellites signals, improving the imaging gain with less computational complexity. In the same year, Nithirochananont et al. [35] explored the passive coherent multistatic SAR using spaceborne illuminators. The k-space support was derived for bistatic SAR images, and the spatial resolution was improved after coherently combination. For this strategy, the improvement of range resolution is based on a multistatic configuration system, and multistatic signals are collected to be processed for the improvement of resolution. This strategy utilizes the rich and diverse resources of GNSS signals, and region of interest is illuminated from different angles. Hence, more targets scattering information from different bistatic SAR imaging geometry can be acquired, and imaging resolution can be enhanced by combining multi-angle SAR images. However, more echo signals need to be recorded and processed in this type of strategy, existing some disadvantages in terms of system cost and signal processing complexity due to the multistatic configuration as a result.
As for the other strategy, many researchers concentrate on the signal processing of passive GNSS-based bistatic SAR to improve the range resolution. In 2015, Ma et al. [36] adopted Galileo E5a-Q and E5b-Q signals together using coherent fusion, and range resolution was improved more than three times by the joint of spectrum at the expense of the maximum detection range. In 2017, Zheng et al. [37] proposed a new range compression method for GNSS-based SAR. The signals of two channels were correlated, and a spectrum equalization method was applied to suppress the side lobes. While the signal noise ratio (SNR) was lost after processing. In 2019, Zeng et al. [38] presented a pre-processing algorithm for Galileo-based bistatic SAR. A dual matched filtering was applied on the joint of band of E5-Q signal, and the range resolution can be improved with side-lobes suppression. In the same year, Zheng et al. [39] proposed a new GPS-based SAR imaging algorithm to improve range resolution. A second-order differentiation (Diff2) operator was applied on the square of range compressed signal with respect to range delay lag, and phase value was reserved by recovery processing. The special property of GPS C/A code correlation makes the Diff2 operator feasible to improve the range performance. For the noise free environment or pointlike target, the range profile of range compression C/A code is standard triangle-like shape, which is available to use the Diff2 operator. However, in a real scenario, this triangle-like shape is not perfect, which could affect the range compression performance, such as additional side lobes. Keeping target side lobes amplitude at a relative low level is quite important for SAR image formation, due to the fact that weak targets may be submerged by high side lobes of strong scatterers, and will be missed in the target detection. Hence, the existence of high side lobes will increase false alarm probability, and affect the weak target image formation. For this second strategy for improving range resolution, it is based on the signal processing of individual passive bistatic SAR, and only one bistatic pair is required compared with the first strategy. Hence, from the view of system cost and signal processing complexity, this type of strategy shows better performance compared with multistatic imaging processing. Although only a certain target scattering information can be obtained in this strategy, it guarantees an improved individual passive bistatic SAR imaging resolution, contributing to a better imaging foundation for multistatic image fusion. So this strategy will be considered to improve the range resolution. In this paper, we will concentrate on the improvement of range compression performance, including a finer range resolution and lower targets side-lobes.
In this paper, an image formation method based on the Diff2 operator is proposed. Comprehensive derivation and analysis of Diff2 operator are investigated, and the reason for generating the side lobes after applying the Diff2 operator on the real experimental data is explored. The part of operator which only affects the range resolution is extracted separately, and the phase recovery is then conducted for azimuth compression. Point targets simulation and real scenario experiment with strong scatters are conducted to validate the proposed imaging algorithm in the aspect of range resolution and targets side-lobes. Furthermore, two types of nearby point targets simulation with identical and nonidentical scattering intensity are conducted and discussed to verify the importance of finer resolution and lower side-lobes performance.
The rest of this paper is organized as follows. The overview of GNSS-based SAR system and general signal processing scheme are presented in Section 2. Section 3 gives the proposed image formation method for improving the range resolution. Simulation and experiment results are shown in Section 4, and discussed in Section 5. Section 6 draws the conclusion.

II. SYSTEM OVERVIEW
This section briefly presents the overview of GNSS-based passive SAR system and general signal processing scheme in the range direction. The spatial geometry is illustrated briefly, and the echo signal model is provided at the same time. As for signal processing, the range compression part is introduced stressfully, because the range resolution improvement method proposed in this paper is applied in this step. The rest of imaging steps, including range cell migration correction and azimuth compression, are similar to the conventional active SAR processing methods.

A. OVERVIEW OF GNSS-BASED PASSIVE SAR
GNSS-based passive SAR utilizes GNSS satellites as the transmitters of opportunity, including GPS, GALILEO, GLONASS and BeiDou. The receiver can be mounted on an aircraft, a land vehicle, or it could be stationary on the ground. The transmitter or moving receiver would provide aperture synthesis, guaranteeing the imaging ability along azimuth direction. In this paper, we assume a ground-based fixed receiver, and synthetic aperture is performed by the moving of GNSS satellite. The system configuration of GNSS-based passive SAR is presented in Fig. 1. The transmitted signals from satellites spread into two types of signals, direct signal and radar signal. The radar signal is the reflected signal from the scenario, and recorded by the radar channel of receiver. The direct signal is used for range compression by cross-correlation with the radar signal, and the radar signal is processed into radar image by imaging algorithm.  The direct channel antenna points towards the satellite directly, which can maximize the direct signal strength, guaranteeing a sufficient SNR in the direct channel. The radar channel antenna points towards ground scenario, recording reflected signals. After collecting the radar signal and the direct signal using reflect channel antenna and direct channel antenna, GNSS tracking processing is applied on the direct signal for the satellite positioning and synchronization, tracking the time delay, phase, and navigation message of the GNSS signal, and multipath propagation effects occurred in the satellite positioning and tracking can be mitigated in this part [40]- [43]. Due to the fact that radar signal and direct signal use the same clocks and local oscillators, so the clock slippage and local oscillator drift is the same among two receiving channels, and the tracking results, including the delay and phase of the direct signal, can be also compensated at the radar signal, and the phase error caused by the receiver clock slippage and local oscillator drift can be eliminated [19]. Then intermediate frequency demodulation is applied on the radar signal, following with SAR twodimensional data formation. For the following step is image formation algorithm. The target position in the reflect channel is determined by the time delay of the round trip from the transmitter to the receiver via the targets in the surrounding terrain. Targets range information can be obtained using SAR range compression technique, and azimuth information can be acquired by the SAR azimuth processing technique using the Doppler history of targets. Finally, focused image will be generated after imaging algorithm processing.
In this paper, we will concentrate on the GNSS-based SAR image formation algorithm, shown in the red dashed rectangle in Fig. 2. The improvement of range resolution and side lobes performance during SAR imaging processing will be investigated.
From the perspective of signal model, the transmitted signals of GNSS can be presented in (1) [19]. Here, GPS signal is adopted in this paper, where C (·) means the C/A code, and D (·) represents the navigation code. f c denotes the carrier frequency, which is at 1575.42MHz for L1 band of GPS. θ 0 means the initial signal phase. The navigation code is used for navigation message, including time information and satellite ephemeris, which will not affect the image formation and will be ignored in the following derivation.
After quadrature demodulation and two-dimensionalization for SAR processing, the direct signal at receiver reference channel can be derived as where tr denotes the fast time, and ta denotes the slow time. R TrRx represents the range between transmitter and receiver, which is time dependent, and can be calculated by denote the location of the transmitter and receiver, respectively. Similarly, the radar signal can be shown as where R TrTar and R RxTar represent the range between transmitter and target, receiver and target, and can be calculated by wherex tar ) denotes the location of the target.

B. OVERVIEW OF GENERAL RANGE PROCESSING FOR GNSS-BASED PASSIVE SAR
Compared with conventional active SAR echo signal, which utilizes chirp signal as the transmitting signal, GNSS-based signal adopts the C/A code, which is a pseudo random noise (PRN) code with 1023 chips. Fig. 3 shows the C/A code and its auto-correlation results. The peak occurs at the delay with zero, and decreases rapidly among other chips. This special property makes the C/A code available for measuring the range accurately. According to the PRN character of the C/A code, range compression can only be achieved by cross-correlation with radar signal, rather than conventional matching filtering technique. The range compressed signal s rc (τ ) can be derived as, where τ d = R TrRx c is the time delay of direct signal, and τ r = (R TrTar + R RxTar ) c is the time delay of radar signal. It can be seen that s rc (t a , τ ) has a basic signal model, We assume that R 0 (t a , τ ) represents the auto-correlation of s 0 (t a , t r ) along the range direction.
Then, s rc (t a , τ ) can be derived as, According to (8), the target range location can be reflected After range compression, the rest of imaging processing is similar to conventional image formation algorithm.
In this paper, the method in the paper [21] is utilized to finish the rest of processing scheme.

III. MODIFIED IMAGE FORMATION ALGORITHM OF GNSS-BASED SAR
In this section, the proposed image formation algorithm for range resolution improvement of GNSS-based SAR is introduced in detail. The range compression algorithm is modified based on the Diff2 method proposed in the current research [39], which is also introduced here for comparison. One-dimensional compression of GPS-based SAR is tested by simulation and experiment using conventional cross-correlation method and Diff2 method, and the emerging problem will be discussed. Then, the Diff2 method is analyzed and separated using mathematical derivation, and a modified method based on the Diff2 is proposed. Finally, the one-dimensional compression is testified by the proposed method.

A. Diff2 METHOD
The Diff2 method is used on the square of range compressed signal [39], and the triangle-like range envelop of compressed signal can be sharpened by using second-order differentiation operator.
In this method, the compressed signal is firstly arithmetical squared [39], shown as, Then, the Diff2 of R 2 with respect to each range compressed pulse delay τ i can be carried out [39], Due to the fact that the Diff2 is appied on the squared range compressed signal, which affects the signal phase used for the azimuth compression. Because the Diff2 method actually only influence the envelop of the range profile, and signal phase is changed as double times of original value, which can be recovered by calculating the half-phase of resulting signal. As for the phase wrapping after applying the square operation, the detected phase will be processed by the azimuth recovery function, and the real phase of range compressed signal can be preserved after recovery processing with a wrapping period of 2π, which will not affect the following azimuth compression. Finally, the range compressed signal after Diff2 method can be expressed as, where H azi is the recovery function of the carrier phase term, and can be expressed as, Theoretically, the assumption of triangle envelop is accurate, and the range profile can be perfectly improved. Fig. 4 shows the performance of Diff2 method using onedimensional simulation of GPS C/A code. The target is located at 1600 range samples, and the envelop is nearly an ideal triangle shape. The range resolution can be effectively improved after Diff2 method, while the side lobes also emerges slightly. It is acceptable in the simulation environment. However, in real scenarios, this envelop is not a standard triangle shape anymore, and differentiation operator is easily influenced by the profile fluctuation, which may lead the operator invalidate.
We conducted a real experiment of GPS-based SAR, which will be performed in the next section in detail. Extracting a pulse from the recorded data, and re-using the Diff2 operator. The comparison result is shown in Fig. 5. In Fig. 5, it can be seen that the range resolution can be improved effectively, while the side lobes spring up around the main lobe, located around 3.1 × 10 4 sample. Hence, in this paper, we will deal with the high unnecessary side lobes generated by the Diff2 method, and derive the modified algorithm in detail.

B. PROPOSED ALGORITHM FOR RANGE COMPRESSION IMPROVEMENT
In order to remove the side lobes caused by the Diff2 operator, the signal model expression of the squared range compressed signal after Diff2 operator is analyzed based on the mathematical derivation. VOLUME 8, 2020 The first-order differentiation of squared signal can be presented as, The first order differentialted signal of the squared range compressed signal is the double production of the range compressed signal and its first order differentiated signal. According to the differential law of function production, the second differentiation is derived as, According to (14), the Diff2 operator has been equally separated into two part, where the first one is dependent of the first order of range compressed signal, and the second one is related to both the range compressed signal and its second order differentiation result.
This equivalent operation in (14) is testified by the onedimensional range compressed signal used in Fig. 5, and the comparison result is illustrated in Fig. 6. According to Fig. 6, the performance of range compressed signal after the equivalent operation is nearly identical to the Diff2 operator. This slight difference results from two reasons. The echo signal data was collected and processed in the digital form, and range signal was sampled by the receiver. Signal performs matrix in the programming implementation, which will lead to the difference. Secondly, differentiation operation was implemented using the MATLAB function 'diff', which is a difference and approximate derivative operation. After differentiation processing, the dimension of range compressed signal would reduce one sample if one time of differentiation was applied. In order to keep the same dimension to adapt to the following processing, we supplemented that missed sample with zero, which would lead to a slight drift. So there exist a little different between Diff2 and equivalent operation. After the equivalent operation, the range main lobes of two lines had the same width, and side lobes were also at the same relative amplitude level, validating that the equivalent transform in (14) is feasible.
Then, in order to analyze the individual function of every component on the range profile, the first part and second part in (14) are both plotted within one figure, together with the sum up result, shown in Fig. 7.   FIGURE 7. Comparison of equivalent operation with its two components using experiment signal.
In Fig. 7, the first part, shown as blue dot dash line, only includes two side lobes, and the location of main lobe performs vacant. These two side lobes will cause the high side lobes in the end. While the second part, shown as red dot line, mainly contributes the main lobe, and the side lobes are lower than that within the black solid line in the meanwhile.
Compared with the width of main lobe in the sum-up result, it was slight broadened near the far side of main lobe, however it can be acceptable upon the obvious improvement for side lobes suppression. In addition, the peak of main lobe in the second part appeared two peaks, rather than one peak as black line, due to the fact that the amplitude envelop of compressed C/A code is triangular, and the sample numbers are limited. It is tolerable during processing one frame only, and the peaks of different frames only slightly fluctuate within the main lobe. So this issue can be addressed by the azimuth accumulation. So if only the second part is adopted in the image formation method, the rang resolution can be improved, and the side lobes can also be suppressed simultaneously. Fig. 8 gives the corresponding comparison result using simulation signal, the similar conclusion can be drawed. A better performance can be acquired in the simulation. The main lobe after using the second part shown as red dash line is identical to the equivalent operation results. The side lobes were suppressed significantly, and almost at zero level.
Hence, in this paper, we utilize the range compressed signal and its second order differentiation result to achieve the range performance improvement, shown as (15), where the main lobe of target can be sharpened, and the side lobes still stay at a lower level, rather than applying the Diff2 operator on the squared range compressed signal. s proposed (t a , τ ) = 2s rc (t a , τ ) · ∂ 2 s rc (t a , τ ) ∂τ i 2 (15) 80338 VOLUME 8, 2020 The first part is the range compressed signal, which can be calcualted by (8). The second part is a differentiation operation applied on the range compressed signal, which is a convolution result, shown in (7) previously. This differentiation can achieved by two ways during the program implementation, listed as follows, • The first way is to apply the differentation on the range compressed signal directly. This way can be easily implemented, while every range line needs to be addtionally processed during the processing.
• The second way is using the relation of convolution and differentation, shown in (16). So the second order differentation of range compressed signal can be achieved by mutiplying the second order differentation of radar signal and original reference signal, or mutiplying the original radar signal and the second order differentation of reference signal.
The latter of the second way will be utilized in our method, due to the fact that the second order differentation of reference signal can be calculated one time only, and then multiplied with the range compressed signal per range lines, shown in (17).
As for the phase recovery for the azimuth compression, which is similar to the derivation of the method in reference [39], the differentiation operation is dependent to range delay τ i only, so it only processes the compressed signal envelop, and will not affect the signal phase after differentiation. The phase recovery function H azi (t a , τ ) actually takes the phase distortion into consider, which is due to the square operation on the range comprressed signal. In the proposed processing method, the multiplication in the second part of (14) will also generates the addition of azimuth phase, which is the sum of phase in the range compressed signal and its second order differentiation. Similarly, this second order differentiation on the range compressed signal will also act on the signal envelop, and the differentiated signal phase is same as orginal phase. So the proposed method also generates the signal phase two times superposition, which is identical to the square operation in Diff2. The rest of imaging algorithm, including direct interference cancellation, range cell migration correction and azimuth compression, can use general image formation method for GNSS-SAR, and we employ the algorithm in reference [21] in this paper. Based on the successfully verified results of one-dimensional simulation and experiment, the two-dimensional imaging will be conducted in the next section to show the improvement of range resolution performance. Fig. 9 provides the flowchart of the proposed algorithm, where the main idea is remarked in the red dashed rectangular box.

IV. SIMULATION AND EXPERIMENT RESULTS
In this section, in order to show the effectiveness of the proposed algorithm in terms of range resolution performance improvement, point targets simulation and stadium experiment will strong scatters of GNSS-based passive SAR are conducted. For comparison, the image formation results of the original conventional processed results will also be provided.

A. POINT TARGETS SIMULATION RESULTS
In the point targets simulation, GNSS-SAR based on GPS L1 signal is adopted as an example. The direct echo signal is VOLUME 8, 2020 constructed by applying a time delay on the C/A code directly, and the radar echo signal is built by calculating the phase delay of the round trip by summing the range from transmitter to receiver via the targets. The simulation parameters are listed in Table 1. The imaging geometry is illustrated in Fig. 10. In the imaging scenario, three point targets are simulated. The coordinates of these targets are (0m,0m), (−300m, −300m), (300m, 300m), and their altitudes are all set at the sea level ground. The receiver is located at (−1000m, 0m, 500m). The transmitter selects one of GPS L1 satellites, centering at 1.0723 × 10 4 km, 2.0639 × 10 4 km, 1.2958 × 10 4 km , which can cooperate as a quasi-monostatic configuration. The dwell time in the processing on the scenario is 100s. The echo signal is processed by the proposed and conventional imaging algorithms, shown in Fig. 11(a) and (b), respectively.
Observing Fig. 11, it can be seen that the range performance is improved significantly using the proposed image formation method, and the range resolution has improved as about 30m, while the traditional one is coarsely about 150m. The azimuth results of these two results is almost the same, and the azimuth focusing performance was not deteriorated after the proposed method, validating that the azimuth phase recovery in our method is effective. Furthermore, additional two simulations of point targets were conducted to verify the importance of the proposed algorithm. The simulation parameters were the same as that listed in Table 1.
The first simulation was conducted to perform the perspective of interest in the range resolution improvement, three nearby targets with same scattering coefficients along the range direction were simulated. The coordinates of three targets were (0m, 0m), (−50m, 15m) and (50m, −15m). The space among these targets was narrower than the traditional range resolution (20m). The imaging results are shown in Fig. 12, where the top sub-image is the result of the proposed algorithm, and the bottom one is the original result. It can be seen that three nearby targets are obscure each other and cannot be distinguished in the original result, while they can be recognized after using the proposed method. Hence, the range resolution improvement of GNSS-based SAR benefits better image targets information, and more details can be revealed.
The second simulation was carried out to show the imaging performance of nearby weak target around strong target. The same parameters and targets position were the same as those in the first compensatory simulation, while the targets scattering intensity were configured in different level. For the center target, the scattering intensity was set as normal, and the other two targets were set as 0.1 compared with the center target intensity. The imaging results are shown in Fig. 13. The top image is the result using the proposed method, and the bottom one is the Diff2 operator result. It can be seen that nearby weak targets can be distinguished clearly after the processing of the proposed method, while they were submerged in the side lobes after using the Diff2 operator. Hence, compared with conventional Diff2 operator, the proposed method can distinguish weak targets around strong targets.

B. REAL SCENARIO EXPERIMENT RESULTS
To further confirm the effectiveness of the proposed algorithm for range performance improvement, a GNSS-based passive SAR real scenario experiment with strong scatters is carried out. The receiver is a dedicated four channels GNSS signal collector. The direct channel antenna used in the experiment is a right-handed circularly polarized antenna. The surveillance channel antenna is a left-handed circularly polarized antenna with 12dB antenna gain. The GPS L5 band signal is adopted. The parameters of the experiment are listed in Table 2. The experiment geometry is given in Fig. 14. The measurement was conducted at 14:10 on 11 May 2019, and the satellite of GPS PRN 30 is selected, which can be formed as a quasi-monostatic configuration. The receiver is fixed at the spectators stand in Beihang University, and the surveillance antenna points to the east. Fig. 15 shows the measuring  equipment, including the reflect channel antenna, direct channel antenna and the receiver. The main imaging area is away from the receiver about 200m-300m. The dwell time in the processing is 20s. Considering the inter-channel interference of the receiver, direct signal is inevitably mixed in the surveillance channel, the CLEAN technique [10], [44], [45] is applied to remove the direction signal interference. Table 3 gives the satellite position and during the whole dwell time.
The imaging results of the proposed algorithm are shown in Fig. 16(a). For comparison, Fig. 16(b) gives the results using conventional imaging algorithm.
The images in Fig. 16 were normalized scaled according to the maximum peak of the images. Compared with the  optical photo in Fig. 14, there are mainly two strong scatters that can be both performed in Fig. 16 (a) and (b), located at range about 250m and 300m, which correspond to the natatorium and the gymnasium, shown in Fig. 17. The imaging performance was improved significantly visually after the processing of the proposed imaging method, especially in the range direction.

V. DISCUSSION
Compared with conventional imaging results, the simulation and experiment results shows the proposed algorithm can achieve a better performance, and the range resolution can be improved effectively. In this section, the processing results of the simulation and the experiment are discussed in detail. The imaging results of the proposed algorithm are also compared with the Diff2 operator method, in order to show the improved imaging ability. Fig. 18 shows the comparison of simulation results by using Diff2 operator method and the proposed method.
Observing Fig. 18, the focusing performance are almost the same, and range resolution in both of them was increased, confirming the feasibility of range compression improvement performance. However, from the view of side lobes in the range direction, it is clearer after using the proposed algorithm in Fig. 18(a), and most of the side lobes were suppressed compared with Fig. 18(b), which exists long trailing along the range direction.
Similarly, the comparison is also provided based on the experiment results, shown in Fig. 19.
In Fig. 19, the side lobes in sub-image (a) was obviously suppressed compared with (b). The whole performance is visually clearer after using the proposed method. Meantime, according to Fig. 16, the range resolution is better than the conventional results. In order to further perform the details of Fig. 19, the range profile at the azimuth center frame was extracted, shown in Fig. 20.
Similar to the conclusion of visual observation, the side lobes were suppressed effectively. The side lobes of the first target at 250m have decreased from 0.27 to 0.08 in terms of normalized amplitude, and that of the second target at 300 are from 0.52 to 0.16, respectively. Here, it is worth noticing that the main lobe of target at range 300m in Fig. 19(a) is slightly deteriorated compared with sub-image (b), due to the fact that we only reserved the second part of Diff2 operator which is useful for both resolution improvement and side lobes suppression. Although the first part of Diff2 operator may achieve the main lobe narrower, it will bring out the side lobes seriously. Overall the imaging performance can be improved after the proposed imaging method.
The computational complexity of the proposed algorithm is also analyzed. The computational load of the range cell migration correction and azimuth compression has been discussed in the reference [21], which needs six range FFTs, four azimuth FFTs and seven complex multiplications, requiring 3N a N r log 2 N r + 2N r N a log 2 N a + 7N a N r complex multiplications. Besides these complex multiplications, the proposed range compression algorithm needs additional three operations, including second order differentiation operation of the range compressed signal, azimuth phase compensation and multiplication operation, requiring 2N a N r + N a N r + N a N r , which sums up to occupy 5N a N r complex multiplications. Hence, the total assumption of the computational load is 3N a N r log 2 N r + 2N r N a log 2 N a + 12N a N r .
The collected echo data was processed offline using the MATLAB software on a ThinkPad laptop with CORE i7-4940MX 3.10GHz processor and 32GB RAM. The time consumed in the proposed imaging algorithm was 527.6s for the 20s dwell time echo data.
A further point is necessary to be discussed here. Similar to the weak point of the Diff2 operator, the image formation algorithm for range compression improvement needs to be considered in a low-SNR environment. The proposed method is also based on the signal multiplication as same as the square operation of range compression signal in Diff2 operator, which will enhance the image contrast. The strong enough targets will be enhanced to a higher degree, while the weak targets are suppressed further. Fig. 21 gives the range profile using traditional cross-correlation method to achieve range compression, there are some weak targets peak at range 100m and 175m about 0.18 of normalized amplitude, which can be detected in this situation. While in Fig. 20(a) and (b), these two targets are very weak, and almost submerged in the noise. This point will be taken into consideration in our future work, and a more feasible image formation method need to be developed for low-SNR scenario.
In order to emphasize the importance of the proposed algorithm, a further analysis with discussion is presented. It can be contributed into three aspects. Firstly, high resolution SAR image formation is very important in the remote sensing field. For the proposed algorithm, better range resolution can be achieved. This improvement has been compared in Fig. 11 and Fig. 16. The range resolution has improved more than two times than the original one, from 21.2m of conventional result to 10.0m for the proposed method. Secondly, high resolution  image formation with low side lobes is another key point in terms of image quality improvement, and weak targets can be detected effectively. Hence, the distinguishability among adjacent targets can be enhanced. In the proposed algorithm, observing the comparison in Fig. 16, Fig. 20 and Fig. 21, lower side lobes can be achieved compared with the conventional algorithm. Take Fig. 20(b) for example, the side lobes of the target at range 300m, peaking at 0.52, are higher than the main lobe of the target at range 250m, peaking at 0.4. If detection threshold was set at a value lower than 0.4 for the successful detection of target at range 250m, the side lobes would be detected at the same time, increasing false alarm probability. However, for Fig. 20(a), the false alarm probability would not increase upon the same detection threshold. Thirdly, the strong side lobes of targets in SAR images will mask the nearby weak targets, and weak targets at the range of strong targets side lobes will be missed, shown in Fig. 12 and Fig. 13, while these weak targets can be recognized using the proposed algorithm.

VI. CONCLUSION
In this paper, a GNSS-based passive SAR image formation has been presented. The reason for causing high side lobes in the range direction using the conventional Diff2 method was analyzed in detail. An improved range processing method based on the Diff2 operator was proposed using mathematical factorization and multiplication of range cross-correlated signal with its second order differenced signal.
The algorithms were verified through simulated point targets and real field GNSS-based SAR experiment. It was shown that after the proposed method processing, the image performance was improved effectively. Selected point targets were compared, and range profiles showed a better range resolution and decreased side lobes, contributing to the better distinguishability of nearby targets and the recognition of weak targets.