An E ﬃ cient Imaging Algorithm for GNSS-R Bi-Static SAR

: Global Navigation Satellite System Reﬂectometry (GNSS-R) based Bi-static Synthetic Aperture Radar (BSAR) is becoming more and more important in remote sensing, given its low power, low mass, low cost, and real-time global coverage capability. Due to its complex conﬁguration, the imaging for GNSS-R BSAR is usually based on the Back-Projection Algorithm (BPA), which is very time consuming. In this paper, an e ﬃ cient and general imaging algorithm for GNSS-R BSAR is presented. A Two Step Range Cell Migration (TSRCM) correction is ﬁrstly applied. The ﬁrst step roughly compensates the RCM and Doppler phase caused by the motion of the transmitter, which simpliﬁes the SAR data into the quasi-mono-static case. The second step removes the residual RCM caused by the motion of the receiver using the modiﬁed frequency scaling algorithm. Then, a cubic phase perturbation operation is introduced to equalize the Doppler frequency modulation rate along the same range cell. Finally, azimuth phase compensation and geometric correction are completed to obtain the focused SAR image. A simulation and experiment are conducted to demonstrate the feasibility of the proposed algorithm, showing that the proposed algorithm is more e ﬃ cient than the BPA, without causing signiﬁcant degradation in imaging quality.


Introduction
Global Navigation Satellite System Reflectometry (GNSS-R)uses the reflected GNSS signal to obtain information about the Earth's surface, such as the speed of ocean winds [1,2], the soil moisture [3][4][5], the biomass of the land [5,6], the ocean altitude [7,8], etc.The concept of GNSS-R was firstly considered by Hall and Cordy in 1988 [9].However, based on their erroneous estimates, they deemed this technique not very promising.In 1993, the concept for ocean altimetry using the GNSS reflected signal was proposed by Martin-Neira [10].From then on, GNSS-R has achieved rapid development, both in the measurement concepts, instruments, and applications [11,12].The first space borne GNSS-R mission, U.K. Disaster Monitoring Constellation (DMC), was launched in 2003 [13].The Low Earth Orbiting U.K. TechDemoSat-1 satellite (TDS-1) was launched in 2014 [1].The Cyclone Global Navigation Satellite System (CYGNSS) aimed to improve the prediction of extreme weather was launched by NASA in 2016 [14].
Firstly, benefiting from the removal of the transmitting device, the GNSS-R system can be made low mass, low power, and low cost.Secondly, we can choose the GNSS satellites with different azimuth and elevation angles to observe the target from different perspectives.Thirdly, it can achieve easy synchronization aided by the GNSS time service [26].Finally, it can achieve nearly real-time global monitoring due to the rich satellite resources.On the other hand, the GNSS-R BSAR also has some disadvantages.Firstly, the power of the GNSS signal is very low around the Earth's surface, and thus, a longer integration time is needed to have a satisfactory Signal-to-Noise Rate (SNR) for the SAR image.Secondly, the instability of the illuminating power will lead to problems with the signal power calibration [27].Another challenge for GNSS-R BSAR lies in efficient imaging.In general, two main difficulties are faced in the imaging of GNSS-R BSAR due to its complex geometry:

•
The complex expression for range history: The range history is due to the motion of both the GNSS satellite and the receiving platform, which makes it extremely difficult to obtain a precise analytical solution to the stationary phase point of the Doppler phase.

•
The translational variant Doppler history: Unlike mono-static SAR, the Doppler history of the echo signal of GNSS-R BSAR is translationally variant since the trajectories of the transmitter and the receiver are non-parallel and their velocities are also different.This means that imaging becomes a two-dimensional spatially varying filtering process.
As a result, commonly used imaging algorithms for GNSS-R BSAR are usually based on the Back-Projection Algorithm (BPA) [15][16][17][18][19][20][21][22][23].Although having the ability of handling different configurations, these algorithms are very time consuming.In order to handle the spatial variation of the echo data, some two step algorithms [24,25] were proposed for GNSS-R BSAR.The general idea of those algorithms is very similar: a bulk compensation for coarse focusing and then a pointwise compensation for the residual RCM and phase error.These algorithms are proven to have better performance than the BPA.However, the second step for residual compensation is still inefficient.In [28], a Nonlinear Chirp Scaling (NCS) algorithm was presented to focus the SAR data collected in the translational variant mode for airborne BSAR.The key step is to equalize the Doppler Frequency Modulation (FM) rate using phase perturbation.This idea could potentially also be extended to the GNSS-R BSAR case.
In this paper, an accurate signal model is firstly established for GNSS-R BSAR.Then, the range history and Doppler phase history of the echo signal are analysed.Based on the analysis, an efficient imaging algorithm for GNSS-R BSAR is proposed.In the first step, the RCM and the Doppler phase caused by the motion of the transmitter are roughly removed through bulk compensation in the azimuth time domain.After this step, the bi-static SAR data are simplified into a quasi-mono-static case, so that precise expression of the echo data in the range-Doppler domain can be derived.In the range-Doppler domain, the RCM caused by the motion of the receiver and the residual RCM caused by the motion of the transmitter are corrected through the modified frequency scaling algorithm.Then, the data are transformed into the azimuth time domain, and phase perturbation is performed to equalize the Doppler FM rate along the same range cell.Finally, the data are transformed into the range-Doppler domain for azimuth phase compensation.In the whole procedure, only complex multiplications and FFT are needed, leading to a very efficient algorithm.
In order to demonstrate the feasibility and efficiency of the proposed imaging algorithm, point target simulation is firstly conducted.The GPS L5 satellite is chosen as the illuminator.The receiver is assumed to be moving along a nonparallel line to the satellite with different velocities.In total, 25 targets are laid out evenly in a 20 km × 20 km square area.The motion error of the receiver is ignored in our simulation.The imaging results of point targets at different locations are very consistent with theoretical analysis, demonstrating the feasibility of the proposed algorithm.
An experiment is also conducted to further confirm the feasibility of the proposed algorithm.Again, the GPS L5 satellite is chosen as the illuminator.The buildings at Beihang University are chosen as the imaging targets.The collected echo data are processed using the traditional BPA and the proposed algorithm, respectively.The imaging results are similar for the two algorithms, but the efficiency of the proposed algorithm is much higher than the traditional BPA.The imaging result is then compared with the optical image, and the strong scattering areas match the high density areas in the radar image well.
The paper is organized as follows: The signal model for GNSS-R BSAR and details of the proposed algorithm are presented in Section 2. Simulation and experimental results are provided in Section 3. A further discussion about the proposed algorithm is presented in Section 4, and conclusions are drawn in Section 5.

Modelling and Analysis
The general geometric configuration of the GNSS-R BSAR is shown in Figure 1.The receiver is mounted on an airplane moving in parallel to the ground along the x direction with a constant velocity V r , and the transmitter is mounted on the GNSS satellite moving along its Kepler orbit.The target is positioned on the ground plane (x-y plane).The bi-static angle θ is defined as the angle between the target-to-satellite baseline and the target-to-receiver baseline.The receiver works in the broadside mode, and the beam centre of the receiver antenna falls at the y axis at t = 0. Here, t represents the azimuth slow time.The beam of the transmitter antenna covers most area of the Earth's surface, and the antenna pattern is considered as constant in the imaging duration.For convenience of expression, in the following part of this paper, the slant range between the target and the satellite is referred to as the transmitter range, and the slant range between the target and the receiver is referred to as the receiver range.The paper is organized as follows: The signal model for GNSS-R BSAR and details of the proposed algorithm are presented in Section 2. Simulation and experimental results are provided in Section 3. A further discussion about the proposed algorithm is presented in Section 4, and conclusions are drawn in Section 5.

Modelling and Analysis
The general geometric configuration of the GNSS-R BSAR is shown in Figure 1.The receiver is mounted on an airplane moving in parallel to the ground along the x direction with a constant velocity Vr, and the transmitter is mounted on the GNSS satellite moving along its Kepler orbit.The target is positioned on the ground plane (x-y plane).The bi-static angle θ is defined as the angle between the target-to-satellite baseline and the target-to-receiver baseline.The receiver works in the broadside mode, and the beam centre of the receiver antenna falls at the y axis at t = 0. Here, t represents the azimuth slow time.The beam of the transmitter antenna covers most area of the Earth's surface, and the antenna pattern is considered as constant in the imaging duration.For convenience of expression, in the following part of this paper, the slant range between the target and the satellite is referred to as the transmitter range, and the slant range between the target and the receiver is referred to as the receiver range.For a typical SAR system, the transmitted signal is Linear Frequency Modulated (LFM).The GNSS system is not designed for SAR imaging application, and the structure of the transmitted signal is different from the LFM wave.For example, the GPS signal has three levels in structure: the bottom level is the carrier wave at the L-band; the middle level is the ranging codes, which are pseudorandom sequences; the upper level is the information code that carries the navigation message.For these three sequences, only the ranging code is used for imaging, and thus, the synchronization process should be performed firstly to cancel the other two sequences [15].
Assume the beam centre of the receiving antenna is focusing at target k at t = t0.Then, the echo data of target k after synchronization and SAR data formation could be expressed as: For a typical SAR system, the transmitted signal is Linear Frequency Modulated (LFM).The GNSS system is not designed for SAR imaging application, and the structure of the transmitted signal is different from the LFM wave.For example, the GPS signal has three levels in structure: the bottom level is the carrier wave at the L-band; the middle level is the ranging codes, which are pseudorandom sequences; the upper level is the information code that carries the navigation message.For these three sequences, only the ranging code is used for imaging, and thus, the synchronization process should be performed firstly to cancel the other two sequences [15].
Assume the beam centre of the receiving antenna is focusing at target k at t = t 0 .Then, the echo data of target k after synchronization and SAR data formation could be expressed as: where a(t) is the antenna pattern of the receiver, C(τ) is the ranging code sequence, τ is the fast time, λ is the wavelength of the carrier, and c is the speed of light.R k (t) is the total slant range, and R k (t) = R trk (t) + R tsk (t).R trk (t) is the receiver range for target k, and R tsk (t) is the transmitter range for target k.
Assume that target k 0 has the same minimum transmitter range as target k.Then: The receiver range has the form of a hyperbola: Using the Equivalent Slant Range Mode (ESRM) [29], the transmitter range can be expressed as: where V 0 is the equivalent velocity and ϕ is the equivalent squint angle.
As can be seen from ( 3), (4), and (5), since the range history of the echo data for GNSS-R BSAR is caused by the motion of both the GNSS satellite and the airplane, the Doppler history of the echo data consists of two square roots.In this case, the stationary phase point does not have an analytical solution [30], and thus, the precise explicit expression of the echo data in the Range-Doppler (R-D) domain cannot be obtained.Although the approximate R-D expression of the echo data for airborne BSAR and some associated imaging algorithms have been developed [31][32][33], the derivation is lengthy and tedious, and the accuracy of the approximation also decreases quickly when the bi-static angle increases.In this paper, a two step method is presented to get the R-D expression of the echo data more accurately for GNSS-R BSAR.It will be proven that the square root, which corresponds to the transmitter range, can be removed firstly by a bulk compensation.Then, the expression of the echo data is simplified into a quasi-mono-static case, which only has one square root.This process is referred to as slant range separation.
The slant range separation process is feasible due to the asymmetry between the transmitter range and the receiver range for GNSS-R BSAR.To show the details, we firstly take a look at the transmitter range.Again, GPS is chosen as an example.
The GPS orbit belongs to the high Earth orbit (HEO).The average altitude is about 20,200 km, and the orbit is very close to a circle.Since the altitude of the satellite is too high, if observed on the Earth's surface, the spatial variation of the Doppler parameters of the GPS satellite is very slow.In fact, when the imaging scene is not too large and the integration time is not too long, the second and higher order components of the spatial variation of the transmitter range are negligible.It will be proven in the following by both theoretical derivation and simulation.
The orbit geometry of the GPS satellite is shown in Figure 2. The Doppler FM rate can be calculated by: Assuming that angle ϕ is constant for targets in the imaging scene, the directional derivative of f r on the ground plane is: where grad(R ts ) is the modulus of the gradient of R ts on the ground plane, which is equal to the cosine of the elevation angle α of the satellite.When ϕ increases, ts decreases, and cosα increases.
When ϕ increases to ϕ max , cosα has the maximum value of one.Therefore: where R ts0 is the height of the satellite relative to the ground when ϕ = 0 and V s is the velocity of the satellite.Then, the maximum range error caused by ignoring the second order component of the spatial variation of the transmitter range is: where T is the integration time and l is the size of the imaging scene.
The corresponding phase error is: The simulation result is shown in Figure 3.The six curves correspond to the maximum phase error as a function of integration time for different sizes of the scene (10 km to 20 km with a linear spacing of 2 km from the bottom to the top).The red straight line corresponds to the π/4 phase error criterion for the GPS L5 carrier.It can be seen that, for an imaging scene no larger than 20 km × 20 km, the phase error caused by the second and higher order components of the spatial variation of the transmitter range is still acceptable even when the integration time reaches 30 s.Thus, if a bulk compensation for the transmitter range using the scene centre as the reference point is firstly applied, the residual transmitter range has only a fixed linear component.Then, the expression of the echo data is simplified into a quasi-mono-static case, which only has one square root, and the stationary phase point can then be calculated analytically.

Remote Sens. 2019, 11, x FOR PEER REVIEW
The corresponding phase error is: The simulation result is shown in Figure 3.The six curves correspond to the maximu error as a function of integration time for different sizes of the scene (10 km to 20 km with spacing of 2 km from the bottom to the top).The red straight line corresponds to the π/4 ph criterion for the GPS L5 carrier.It can be seen that, for an imaging scene no larger than 20 km × 20 km, the phase error caused by the second and higher order components of th variation of the transmitter range is still acceptable even when the integration time reaches 30 seconds.Thus, if a bulk compensation for the transmitter range using the scene cent reference point is firstly applied, the residual transmitter range has only a fixed linear com Then, the expression of the echo data is simplified into a quasi-mono-static case, which only square root, and the stationary phase point can then be calculated analytically.where Rtsc(t) is the transmitter range for target c, which is positioned at the centre of the scen Expanding   −   at t = 0 and keeping only those up to the first order, we ha of the transmitter range is still acceptable even when the integration time reaches ds.Thus, if a bulk compensation for the transmitter range using the scene centre as the point is firstly applied, the residual transmitter range has only a fixed linear component.e expression of the echo data is simplified into a quasi-mono-static case, which only has one ot, and the stationary phase point can then be calculated analytically.r the bulk compensation for transmitter range, the echo signal for target k0 becomes: 11) tsc(t) is the transmitter range for target c, which is positioned at the centre of the scene.anding   −   at t = 0 and keeping only those up to the first order, we have: tk is the position shift of the signal in the echo window.The Doppler frequency fd and its e can be expressed as a function of the parameters of the satellite and the target.
⃗ is the unit directional vector of the x axis,  ⃗ is the unit directional vector of the y axis, and e unit directional vector from the target to the satellite.stituting ( 12) into (11) leads to: Phase error/rad After the bulk compensation for transmitter range, the echo signal for target k 0 becomes: where R tsc (t) is the transmitter range for target c, which is positioned at the centre of the scene.Expanding R tsk (t) − R tsc (t) at t = 0 and keeping only those up to the first order, we have: where: R shiftk is the position shift of the signal in the echo window.The Doppler frequency f d and its derivative can be expressed as a function of the parameters of the satellite and the target.
where → i x is the unit directional vector of the x axis, → i y is the unit directional vector of the y axis, and → i rts is the unit directional vector from the target to the satellite.
Substituting ( 12) into (11) leads to: where: The last approximation holds when the scene size and integration time are not too large.
where k a is the Doppler FM rate for target k caused by the motion of the receiver.The phase ϕ 1 is omitted because it has no impact on the derivation.For mono-static SAR, the echo data for targets with the same minimum slant range have the same Doppler history and RCM curve.This means the echo data are translationally invariant.Thus, the imaging algorithms in the azimuth frequency domain are usually more efficient.Due to the non-parallel trajectory and different velocity of the transmitter and the receiver, the echo data for GNSS-R BSAR are translationally variant.This property can be observed by (17) in two aspects.

1.
The spatially varying Doppler centroid: It is formed by the linear component of the residual transmitter range and represented by the spatially varying delay time t dk .This means the equivalent squint angle is also spatially varying.

2.
The translational variant Doppler FM rate: It is formed by the constant component of the residual transmitter range and represented by the shifting factor R shiftk .The shifting factor R shiftk indicates that the echo data for targets with the same minimum receiver range (thus the same Doppler FM rate) will not appear in the same range cell.In other words, the signals that appear in the same range cell have a different Doppler FM rate.
To illustrate the range cell shift effect of the shifting factor, two spaces are defined for GNSS-R BSAR, as shown in Figure 4.

•
The scene space in which the echo data are received.

•
The echo space in which the echo data are stored and processed.
Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 22 where: The last approximation holds when the scene size and integration time are not too large.
where ka is the Doppler FM rate for target k caused by the motion of the receiver.The phase φ1 is omitted because it has no impact on the derivation.For mono-static SAR, the echo data for targets with the same minimum slant range have the same Doppler history and RCM curve.This means the echo data are translationally invariant.Thus, the imaging algorithms in the azimuth frequency domain are usually more efficient.Due to the non-parallel trajectory and different velocity of the transmitter and the receiver, the echo data for GNSS-R BSAR are translationally variant.This property can be observed by (17)  To illustrate the range cell shift effect of the shifting factor, two spaces are defined for GNSS-R BSAR, as shown in Figure 4.

•
The scene space in which the echo data are received.

•
The echo space in which the echo data are stored and processed.As can be seen from the figure, due to the effect of the shifting factor Rshiftk, for targets with the same minimum receiver range (thus the same Doppler FM rates since they are caused only by the motion of the receiver after bulk compensation for transmitter range) in the scene space, their echo data will not appear in the same range cell in the echo space.In other words, for the echo data in the same range cell, their Doppler FM rates are different.As an example, target k has the same minimum As can be seen from the figure, due to the effect of the shifting factor R shiftk , for targets with the same minimum receiver range (thus the same Doppler FM rates since they are caused only by the motion of the receiver after bulk compensation for transmitter range) in the scene space, their echo data will not appear in the same range cell in the echo space.In other words, for the echo data in the same range cell, their Doppler FM rates are different.As an example, target k has the same minimum receiver range as target k 0 and thus has the same Doppler FM rate as target k 0 .However, the echo data of target k appear in the same range cell as target k 1 rather than target k 0 in the echo space, since R k In order to focus the echo data efficiently in the R-D domain, the Doppler FM rate in the same range cell needs to be equalized.Here, the term "equalize" means making the Doppler FM rates equal.Firstly, we consider the shifting curve R shiftk .
The first term of R shiftk is caused by the difference in transmitter range for different targets at t = 0.This term is positively related to the bi-static angle θ.Since θ varies widely in different receiving geometry, this term also varies widely and in general is also the dominating term of R shiftk .This term is a linear function of the azimuth time.The second term is caused by the variation of the Doppler frequency along the direction of the receiver velocity, and this is a quadratic function of the azimuth time.The third term is caused by the variation of the Doppler frequency along the cross direction of the receiver velocity, and it is also a linear function of the azimuth time, but the slope varies with y 0 .
The latter two terms are a function of x 0 and y 0 , which are negligible when the imaging scene is not too large.Ignoring the second term, we have: which is a linear function of the azimuth time.
Correspondingly, for targets whose echo data appear in the same range cell, the minimum receiver range is: where: Then, the Doppler FM rate along the same range cell is: It can be seen that the Doppler FM rate along the same range cell is approximately a linear function of the azimuth time.Therefore, a cubic perturbation function is sufficient to equalize the Doppler FM rate along the range cell.

Imaging Methods
Based on the above analysis, an efficient and general imaging algorithm is proposed for GNSS-R BSAR.The block diagram of the proposed algorithm is shown in Figure 5. Firstly, two receiving channels are used to collect the reflected signal and the direct signal, respectively.The reflected channel signal is used for imaging, whereas the direct channel signal is used for positioning and synchronization.This step is introduced to track the time delay, phase, and navigation message of the GNSS signal and eliminate the phase error caused by clock slippage and local oscillator drift of the receiver [15].
Then, the signal is truncated to form the SAR data.The echo data of target k after synchronization and SAR data formation can be expressed as (11).

Bulk Compensation and Range Compression
Fast Fourier Transform (FFT) along the range direction is performed, and bulk compensation for transmitter range along with range compression is accomplished in the range frequency domain.The expression of the bulk compensation and range compression filter function is given by:  =  *   2   (26) where  *  is the complex conjugate of the Fourier transform of C(τ).Then, the data are transformed into the time domain by performing the range IFFT, and the signal becomes: where P(τ) is the auto-correlation function of the GNSS ranging code.The data are then transformed into the range-Doppler domain.

Residual RCMC
Due to the different equivalent squint angle caused by the linear component of the residual transmitter range, the Range Cell Migration (RCM) of the data in the same range cell has different slope in the azimuth time domain.However, they should still satisfy the same RCM curve in the range-Doppler domain, just corresponding to different intervals, and thus can be corrected Firstly, two receiving channels are used to collect the reflected signal and the direct signal, respectively.The reflected channel signal is used for imaging, whereas the direct channel signal is used for positioning and synchronization.This step is introduced to track the time delay, phase, and navigation message of the GNSS signal and eliminate the phase error caused by clock slippage and local oscillator drift of the receiver [15].
Then, the signal is truncated to form the SAR data.The echo data of target k after synchronization and SAR data formation can be expressed as (11).

Bulk Compensation and Range Compression
Fast Fourier Transform (FFT) along the range direction is performed, and bulk compensation for transmitter range along with range compression is accomplished in the range frequency domain.The expression of the bulk compensation and range compression filter function is given by: where F * ( f τ ) is the complex conjugate of the Fourier transform of C(τ).Then, the data are transformed into the time domain by performing the range IFFT, and the signal becomes: where P(τ) is the auto-correlation function of the GNSS ranging code.The data are then transformed into the range-Doppler domain.

Residual RCMC
Due to the different equivalent squint angle caused by the linear component of the residual transmitter range, the Range Cell Migration (RCM) of the data in the same range cell has different slope in the azimuth time domain.However, they should still satisfy the same RCM curve in the range-Doppler domain, just corresponding to different intervals, and thus can be corrected uniformly according to the reference point, which is chosen to have a zero equivalent squint angle and thus a zero delay time t d .
In traditional SAR applications, the chirp scaling algorithm is often used to correct the range cell variant RCM of the SAR data.However, in GNSS-R BSAR, the transmitted signal is not a chirp signal, and thus, the chirp scaling algorithm is no longer suitable.The Frequency Scaling (FS) algorithm was introduced in [34] to correct the RCM of the de-ramped signal for spotlight SAR.In this paper, it will be proven that after some changes, this algorithm is also capable of handling the RCM Correction (RCMC) problem for GNSS-R BSAR:

•
The scaling operation is not applied to the range frequency axis, but to the range time axis, and thus in some sense, it is more appropriate to name the modified algorithm the range scaling algorithm rather than the frequency scaling algorithm.

•
The scaling of the range time needs four multiplications rather than three as in the frequency scaling algorithm.

•
The scaling factors should all be adapted to the new application.
The detail of the residual RCMC process is illustrated in the following.Firstly, the RCM caused by the motion of the receiver and the residual RCM caused by the motion of the transmitter are equalized by the so-called range scaling algorithm.This step is a combination of four sequential processes: frequency shift, range scaling, inverse frequency shift, and inverse range scaling.After these four operations, the RCM of the signal at different range cells is equalized, and then, a bulk RCMC is applied to finish correcting the RCM.
Considering the reference point for target k 0 (note that the reference point is not always located at the y axis, and target k 0 has a non-zero equivalent squint angle), the echo data of the reference point after bulk compensation for transmitter range have the following form: The RCM in the azimuth time domain is given by: Transformed into the range-Doppler domain, the RCM becomes: where: R trk0 = R trk0 (0) Thus, the RCMC can be accomplished by two main steps: scaling the range time axis by β times, followed by a bulk RCMC to correct the equalized RCM αβR trc .The four scaling filter functions corresponding to the first step are: Here, the subscripts in the four filter functions represent frequency shift, range scaling, inverse frequency shift, and inverse range scaling, respectively.The first multiplication is performed to shift the range frequency of the echo data according to their RCM.After this, the base-band range signal is modulated by a different frequency according to their RCM.The 2D frequency spectrum has a similar curve as the RCM along the azimuth direction, and the range scaling multiplication is then performed.The factor k in the frequency shift filter function H fs is not fixed, but it should make sure that the spectrum shift does not surpass the sampling rate.The purpose of the next three filter functions is similar to their counterparts in the frequency scaling algorithm.For a detailed description of their functions, please refer to [34].The final bulk RCMC filter function is given by: This filter function can be incorporated into the inverse frequency shift filter function in practice, given by: Note that in our analysis, we removed the range cell shift of the targets from the reference point.This will cause a small error in the RCM correction for those targets.However, since the range cell shift is so small in terms of RCM spatial variation, this error is actually negligibly small.

Azimuth Phase Perturbation
After residual RCMC, the RCM of the echo data is corrected.The data are then transformed into the azimuth time domain.The expression of the data after residual RCMC for target k becomes: Since β is very close to one, the broadening of the range profile is negligible, and then: As shown in our earlier analysis, the cubic phase perturbation function is sufficient to equalize the Doppler FM rate along the range cell.Considering the expression of ∆k, the expression of the perturbation function should be: In addition to equalizing the Doppler FM rate along the range cell, the perturbation filter function also has some side effects.The effect of the components with different orders in the perturbation filter function is analysed below.

•
The constant component πat 0 3 : It changes the phase of the signal according to the position of the target, which can be ignored when only the amplitude of the image is concerned.Here: • The first order component 3πat 0 2 (t − t 0 ): It adds a small spatially varying Doppler shift to the signal.This shift is a quadratic function of the azimuth time t, which can be incorporated into the Doppler shift caused by the residual transmitter range, i.e., (41)

•
The second order compound 3πat 0 (t − t 0 ) 2 : It equalizes the Doppler FM rate along the range cell.

•
The third order compound πa(t − t 0 ) 3 : It is a cubic phase modulation, which is the same for all targets.It is far smaller than the phase modulation caused by the receiver motion, which can be ignored during the derivation of the azimuth stationary phase point of the signal.

Azimuth Compensation
After applying the phase perturbation, the second and higher order components of the echo data along the same range cell are equalized, and thus, those data can be focused by the same azimuth compensation filter function.Considering the reference target k 0 , the expression of the echo data of target k 0 after phase perturbation is: where: The method of the stationary phase is applied to calculate the expression of the echo data for the reference target k 0 in the range-Doppler domain.Ignoring the cubic phase modulation caused by the phase perturbation during the derivation of the stationary phase point, we have: The azimuth compensation filter function is chosen as: The azimuth compensation is accomplished by multiplication of the echo signal with the azimuth compensation filter function.Then, the data are transformed into the time domain by the inverse FFT to obtain the focused SAR image.For target k, it is: where:

Simulation Results
Simulations were conducted to demonstrate the performance of the proposed algorithm.The geometric relations of the satellite and the receiver are shown in Table 1.Some parameters used in the simulation are listed in Table 2.An array of 25 targets were considered, which were laid out on a square area of 20 km × 20 km, as shown in Figure 6a, and the imaging results using the proposed algorithm are shown in Figure 6b.The bi-static angle θ was about 25 degrees, and the angle between the target velocity and the transmitter velocity was about 26 degrees.A rectangular antenna pattern was assumed in both the azimuth and range directions.The GPS L5 satellite was chosen as the illuminator, and the ranging code sequence of the GPS L5 signal and its auto-correlation sequence are shown in Figure 7.The imaging results of five chosen targets (Targets 1, 7, 13, 19, and 25) are presented in Figure 8.As can be seen, Targets 7, 13, and 19 were well focused.The range profile and the azimuth profile were well consistent with theoretical ones.For Targets 1 and 25, slight degradation can be observed in the azimuth profile.This is also consistent with our analysis, since Targets 1 and 25 were positioned at the edge of the scene, where the quadratic part of the range cell shift reached its maximum value (about 184 m in the simulation).Thus, Targets 1 and 25 suffered from the most severe Doppler FM rate mismatch, which resulted in distortion in the azimuth profile.The evaluation results are listed in Table 3.For Targets 7, 13, and 19, the degradation in the Peak Side-Lobe Ratio (PSLR) and Integrated Side-Lobe Ratio (ISLR) was in fact negligible.For Targets 1 and 25, the degradation in azimuth PSLR was 0.29 dB and 0.36 dB, respectively, and the degradation in azimuth ISLR was 0.40 dB and 0.31 dB, respectively.The GPS L5 satellite was chosen as the illuminator, and the ranging code sequence of the GPS L5 signal and its auto-correlation sequence are shown in Figure 7.The GPS L5 satellite was chosen as the illuminator, and the ranging code sequence of the GPS L5 signal and its auto-correlation sequence are shown in Figure 7.The imaging results of five chosen targets (Targets 1, 7, 13, 19, and 25) are presented in Figure 8.As can be seen, Targets 7, 13, and 19 were well focused.The range profile and the azimuth profile were well consistent with theoretical ones.For Targets 1 and 25, slight degradation can be observed in the azimuth profile.This is also consistent with our analysis, since Targets 1 and 25 were positioned at the edge of the scene, where the quadratic part of the range cell shift reached its maximum value (about 184 m in the simulation).Thus, Targets 1 and 25 suffered from the most Amplitude/dB The imaging results of five chosen targets (Targets 1, 7, 13, 19, and 25) are presented in Figure 8.As can be seen, Targets 7, 13, and 19 were well focused.The range profile and the azimuth profile were well consistent with theoretical ones.For Targets 1 and 25, slight degradation can be observed in the azimuth profile.This is also consistent with our analysis, since Targets 1 and 25 were positioned at the edge of the scene, where the quadratic part of the range cell shift reached its maximum value (about 184 m in the simulation).Thus, Targets 1 and 25 suffered from the most severe Doppler FM rate mismatch, which resulted in distortion in the azimuth profile.The evaluation results are listed in Table 3.For Targets 7, 13, and 19, the degradation in the Peak Side-Lobe Ratio (PSLR) and Integrated Side-Lobe Ratio (ISLR) was in fact negligible.For Targets 1 and 25, the degradation in azimuth PSLR was 0.29 dB and 0.36 dB, respectively, and the degradation in azimuth ISLR was 0.40 dB and 0.31 dB, respectively.An azimuth variant rotation in the range side-lobe can be observed in the 2D contour of the imaging results.This was caused by the linear component of the spatially varying residual transmitter range.It can also be considered as the system having a spatially varying equivalent squint angle.

Remote
The azimuth resolution for GNSS-R BSAR was mainly determined by the receiver due to the asymmetrical geometry.For the range resolution, however, the analysis was quite complex [35].In addition to the signal bandwidth and the grazing angle of the transmitter and the receiver, the bi-static angle also had a great impact on the range resolution.Generally speaking, the range resolution became worse when the bi-static angle increased.When the transmitter and the receiver were located at different sides of the imaging scene, the range resolution became very poor.Thus, this receiving geometry should be avoided in practice.

Experimental Results
To further confirm the feasibility of the proposed algorithm, a GNSS based BSAR experiment was conducted.Some experimental setups are shown in Figure 9, with some of the parameters listed in Table 4.The receiver was a dedicated GNSS signal collector with four receiving channels.The direct channel antenna was a right handed circularly polarized antenna, while the reflect channel signal was a left handed circularly polarized antenna with a 12 dB gain.An azimuth variant rotation in the range side-lobe can be observed in the 2D contour of the imaging results.This was caused by the linear component of the spatially varying residual transmitter range.It can also be considered as the system having a spatially varying equivalent squint angle.
The azimuth resolution for GNSS-R BSAR was mainly determined by the receiver due to the asymmetrical geometry.For the range resolution, however, the analysis was quite complex [35].In addition to the signal bandwidth and the grazing angle of the transmitter and the receiver, the bi-static angle also had a great impact on the range resolution.Generally speaking, the range resolution became worse when the bi-static angle increased.When the transmitter and the receiver were located at different sides of the imaging scene, the range resolution became very poor.Thus, this receiving geometry should be avoided in practice.

Experimental Results
To further confirm the feasibility of the proposed algorithm, a GNSS based BSAR experiment was conducted.Some experimental setups are shown in Figure 9, with some of the parameters listed in Table 4.The receiver was a dedicated GNSS signal collector with four receiving channels.The direct channel antenna was a right handed circularly polarized antenna, while the reflect channel signal was a left handed circularly polarized antenna with a 12 dB gain.The experiment was conducted at 14:10 on 11 May 2019.The buildings at Beihang University were chosen as the targets, and the receiver was positioned at the rostrum of the football field, which was three meters above the ground.The reflect channel antenna was overlooking the area to the east of the rostrum.

Receiver
The GPS satellite PRN30 was chosen as the optimal radio source, positioned at the west of the receiver during the experiment, and was moving southward.The position and velocity of the chosen satellite at the start, middle, and end of measurement calculated by the precise orbit are given in Table 5.The tracked Doppler frequency of the direct channel signal is shown in Figure 10.The Doppler spectrum of the reflect channel signal is shown in Figure 11.The collected echo was processed using the proposed algorithm, and the radar image with 300 s integration time was obtained successfully.The resultant radar image and the optical image of the experimental area are shown in Figure 12 together for comparison.Several strong scattering areas could be found in the image.They were the wire meshes on the east of the football filed, the wire meshes on the east of the volleyball court, the edge of the natatorium, the edge of the gymnasium, the spire of the gymnasium, and the edge of the teaching building (marked as 1-6 in the figure).It can be seen that the positions of those strong scattering areas were well consistent in the two images, which proved the feasibility of the proposed algorithm.151 Hz The GPS satellite PRN30 was chosen as the optimal radio source, positioned at the wes receiver during the experiment, and was moving southward.The position and velocity chosen satellite at the start, middle, and end of measurement calculated by the precise or given in Table 5.The tracked Doppler frequency of the direct channel signal is shown in Fig The Doppler spectrum of the reflect channel signal is shown in Figure 11.The collected ech processed using the proposed algorithm, and the radar image with 300 s integration tim obtained successfully.The resultant radar image and the optical image of the experimental a shown in Figure 12 together for comparison.Several strong scattering areas could be found image.They were the wire meshes on the east of the football filed, the wire meshes on the the volleyball court, the edge of the natatorium, the edge of the gymnasium, the spire gymnasium, and the edge of the teaching building (marked as 1-6 in the figure).It can be se the positions of those strong scattering areas were well consistent in the two images, which the feasibility of the proposed algorithm.

Discussion
In this section, the computational complexity of the proposed algorithm is analysed and compared with the BPA.
In the proposed algorithm, six range FTs, four azimuth FTs, and seven complex multiplications were performed.The computational load can be derived as 3     2   + 2     2   + 7    , where Na and Nr are the number of range and azimuth samples, respectively.The BPA entails two range FTs and one complex multiplication in the range compression step, and its computational load is      2   +     .The computational load in the back-projection step can be expressed as         , where Nx and Ny are the numbers of samples of the imaging scene, and Nkel is the length of the interpolation kernel.Thus, the total computational load for the BPA can be represented by      2   +     +         .Clearly, when the size of the imaging scene increases, the computational load will increase rapidly for the BPA.
To verify our analysis, the collected data were processed by the proposed algorithm and the BPA, respectively, carried out in MATLAB (MathWorks, USA) on a computer with Xeon E3-1535M

Discussion
In this section, the computational complexity of the proposed algorithm is analysed and compared with the BPA.
In the proposed algorithm, six range FTs, four azimuth FTs, and seven complex multiplications were performed.The computational load can be derived as 3N a N r log 2 N r + 2N r N a log 2 N a + 7N a N r , where N a and N r are the number of range and azimuth samples, respectively.The BPA entails two range FTs and one complex multiplication in the range compression step, and its computational load is N a N r log 2 N r + N a N r .The computational load in the back-projection step can be expressed as N kel N x N y N a , where N x and N y are the numbers of samples of the imaging scene, and N kel is the length of the interpolation kernel.Thus, the total computational load for the BPA can be represented by N a N r log 2 N r + N a N r + N kel N x N y N a .Clearly, when the size of the imaging scene increases, the computational load will increase rapidly for the BPA.
To verify our analysis, the collected data were processed by the proposed algorithm and the BPA, respectively, carried out in MATLAB (MathWorks, USA) on a computer with Xeon E3-1535M 2.9 GHz processors and 64 GB RAM.The imaging results are shown in Figure 13.It can be seen that the two images were very similar to each other.The parameters used in the processing are listed in Table 6.The time consumed by the BPA was about 55.6 times more than the proposed algorithm by our analysis.

Na
Nr Nx Ny N kel 3e5 1024 500 500 8 In the experiment, the time consumed in the imaging step using the proposed algorithm was 64 s, while that of the BPA was 6624 s.The ratio between the time consumed by the two algorithms was consistent with the theoretical value.When the size of the scene is larger, this ratio would increase, and the efficiency of the proposed algorithm would be even better than the BPA.
Further, the profiles of the radar image of the gymnasium using these two algorithms are presented in Figure 14a and Figure 14b, respectively, where the profiles acquired using both algorithms matched each other well, with almost the same high imaging accuracy.In the experiment, the time consumed in the imaging step using the proposed algorithm was 64 s, while that of the BPA was 6624 s.The ratio between the time consumed by the two algorithms was consistent with the theoretical value.When the size of the scene is larger, this ratio would increase, and the efficiency of the proposed algorithm would be even better than the BPA.
Further, the profiles of the radar image of the gymnasium using these two algorithms are presented in Figure 14a,b, respectively, where the profiles acquired using both algorithms matched each other well, with almost the same high imaging accuracy.
was consistent with the theoretical value.When the size of the scene is larger, this ratio would increase, and the efficiency of the proposed algorithm would be even better than the BPA.
Further, the profiles of the radar image of the gymnasium using these two algorithms are presented in Figure 14a and Figure 14b, respectively, where the profiles acquired using both algorithms matched each other well, with almost the same high imaging accuracy.

Conclusions
An efficient imaging algorithm was proposed for GNSS-R bi-static SAR.It consisted of a two step RCMC process to acquire an accurate analytical solution to the stationary phase point and an azimuth phase perturbation step to equalize the Doppler FM rate.As demonstrated by point target simulation results, the proposed algorithm performed well in terms of both imaging quality and efficiency.An experiment using the GPS L5 satellite as the illuminator and the buildings at Beihang University as the targets was also carried out, and the resultant image using the proposed algorithm was compared with the optical image and the radar image obtained by the BPA.As shown, the proposed algorithm performed far better than the traditional BPA in terms of efficiency, with no significant degradation of imaging quality.

Figure 2 .
Figure 2. Orbit model of the GPS satellite.Figure 3. Maximum phase error as a func integration time T and scene size l.

Figure 3 .
Figure 2. Orbit model of the GPS satellite.Figure 3. Maximum phase error as a func integration time T and scene size l.

Figure 2 .
Figure 2. Orbit model of the GPS satellite.

ure 2 .
Orbit model of the GPS satellite.

Figure 3 .
Maximum phase error as a function of integration time T and scene size l.

Figure 3 .
Figure 3. Maximum phase error as a function of integration time T and scene size l.
in two aspects.1.The spatially varying Doppler centroid: It is formed by the linear component of the residual transmitter range and represented by the spatially varying delay time tdk.This means the equivalent squint angle is also spatially varying.2. The translational variant Doppler FM rate: It is formed by the constant component of the residual transmitter range and represented by the shifting factor Rshiftk.The shifting factor Rshiftk indicates that the echo data for targets with the same minimum receiver range (thus the same Doppler FM rate) will not appear in the same range cell.In other words, the signals that appear in the same range cell have a different Doppler FM rate.

Figure 4 .
Figure 4. Illustration of the range cell shift effect.(a) is the scene space, and (b) is the signal space.Five targets and their echoes are shown in the figure; the curvature of the echo is proportional to the Doppler Frequency Modulation (FM) rate of the target.

Figure 4 .
Figure 4. Illustration of the range cell shift effect.(a) is the scene space, and (b) is the signal space.Five targets and their echoes are shown in the figure; the curvature of the echo is proportional to the Doppler Frequency Modulation (FM) rate of the target.

Figure 5 .
Figure 5.A block diagram of the proposed algorithm.

Figure 5 .
Figure 5.A block diagram of the proposed algorithm.
An array of 25 targets were considered, which were laid out on a square area of 20 km × 20 km, as shown in Figure6a, and the imaging results using the proposed algorithm are shown in Figure6b.The bi-static angle θ was about 25 degrees, and the angle between the target velocity and the transmitter velocity was about 26 degrees.A rectangular antenna pattern was assumed in both the azimuth and range directions.

Figure 6 .
Figure 6.Layout of the point targets in the simulation and the imaging results.(a) is the layout of the point targets, and (b) is the imaging results.

Figure 7 .
Figure 7.The ranging code sequence of PRN30 and its auto-correlation result.(a) is the ranging code sequence, and (b) is the auto-correlation result.

Figure 6 .
Figure 6.Layout of the point targets in the simulation and the imaging results.(a) is the layout of the point targets, and (b) is the imaging results.
Sens. 2019, 11,  x FOR PEER REVIEW 12 An array of 25 targets were considered, which were laid out on a square area of 20 km × 20 km, as shown in Figure6a, and the imaging results using the proposed algorithm are shown in Figure6b.The bi-static angle θ was about 25 degrees, and the angle between the target velocity and the transmitter velocity was about 26 degrees.A rectangular antenna pattern was assumed in both the azimuth and range directions.

Figure 6 .
Figure 6.Layout of the point targets in the simulation and the imaging results.(a) is the layout of the point targets, and (b) is the imaging results.

Figure 7 .
Figure 7.The ranging code sequence of PRN30 and its auto-correlation result.(a) is the ranging code sequence, and (b) is the auto-correlation result.

Figure 7 .
Figure 7.The ranging code sequence of PRN30 and its auto-correlation result.(a) is the ranging code sequence, and (b) is the auto-correlation result.

Figure 8 .
Figure 8. Evaluated results for the chosen targets: each subfigure has four parts, including the azimuth profile (upper left), range profile (left bottom), 3D contour (upper right), and 2D contour (right bottom) of the chosen target.(a-e) are the evaluated results for target 1, 7, 13, 19, 25, respectively.

Figure 8 .
Figure 8. Evaluated results for the chosen targets: each subfigure has four parts, including the azimuth profile (upper left), range profile (left bottom), 3D contour (upper right), and 2D contour (right bottom) of the chosen target.(a-e) are the evaluated results for target 1, 7, 13, 19, 25, respectively.
Optical image of the experimental area.The transmitter is positioned at the west of the receiver, while the imaging area is positioned at the east of the receiver.(b)Reflect channel antenna and the targets.

Figure 9 .
Figure 9.The experimental setups.(a) is the optical image of the experimental area.(b) is the reflect channel antenna and the targets.(c) is the direct channel antenna.(d) is the receiver.

Figure 9 .
Figure 9.The experimental setups.(a) is the optical image of the experimental area.(b) is the reflect channel antenna and the targets.(c) is the direct channel antenna.(d) is the receiver.

Figure 10 .
Figure 10.The tracked carrier frequency of the direct channel signal.

Figure 11 .
Figure 11.The Doppler spectrum of th channel signal for 300 seconds.

Figure 10 .
Figure 10.The tracked carrier frequency of the direct channel signal.

Figure 11 .
Figure 11.The Doppler spectrum of the reflect channel signal for 300 seconds.Figure 11.The Doppler spectrum of the reflect channel signal for 300 s.

Figure 11 .
Figure 11.The Doppler spectrum of the reflect channel signal for 300 seconds.Figure 11.The Doppler spectrum of the reflect channel signal for 300 s.Remote Sens. 2019, 11, x FOR PEER REVIEW 18 of 22

Figure 12 .
Figure 12.Comparison between the optical image and the resultant radar image.(a) is the optical image, and (b) is the radar image.

Figure 12 .
Figure 12.Comparison between the optical image and the resultant radar image.(a) is the optical image, and (b) is the radar image.
Remote Sens. 2019, 11, x FOR PEER REVIEW 19 of 22 (a) Imaging result using the proposed algorithm.(b) Imaging results using the BPA.

Figure 23 .
Figure 23.Imaging results using the proposed algorithm and the Back Projection Algorithm (BPA).The size of the imaging area is 100 m × 500 m, and the integration time is 300 s.(a) is the imaging results acquired using the proposed algorithm, and (b) is the imaging results acquired using the BPA.

Figure 13 .
Figure 13.Imaging results using the proposed algorithm and the Back Projection Algorithm (BPA).The size of the imaging area is 100 m × 500 m, and the integration time is 300 s.(a) is the imaging results acquired using the proposed algorithm, and (b) is the imaging results acquired using the BPA.

Figure 34 Figure 14 .
Figure34 The azimuth profile of the resultant radar image of the gymnasium.The blue line represents the result using the proposed algorithm, while the red line represents the result by the BPA.(a) is the azimuth profile, and (b) is the range profile.

Table 1 .
The geometric relations of the simulation.

Table 2 .
Parameters used in the simulation.

Table 3 .
Evaluation results for the chosen targets.ISLR, Integrated Side-Lobe Ratio.

Table 4 .
Parameters used in the experiment.

Table 5 .
The position and velocity of the satellite.

Table 4 .
Parameters used in the experiment.

Table 5 .
The position and velocity of the satellite.

Table 6 .
Parameters used in the processing.

Table 6 .
Parameters used in the processing.