A Uniﬁed Algorithm for the Sliding Spotlight and TOPS Modes Data Processing in Bistatic Conﬁguration of the Geostationary Transmitter with LEO Receivers

: This paper deals with the imaging problem for sliding spotlight (SS) and terrain observation by progressive scan (TOPS) modes in bistatic conﬁguration of the geostationary (GEO) transmitter with a low earth orbit satellite (LEO) receiver, named GTLR-BiSAR system. A uniﬁed imaging algorithm is proposed to process the GTLR-BiSAR data acquired in SS or TOPS modes. Our main contributions include four aspects. Firstly, the imaging geometry of this novel conﬁguration is described in detail. Furthermore, the GTLR-BiSAR signal expressions were deduced in both time and frequency domains. These signal expressions provide great support for the design of processing the algorithm theoretically. Secondly, we present a uniﬁed deramping-based technique according to the special geometry of GTLR-BiSAR to overcome the azimuth spectrum aliasing phenomenon, which typically affects SS and TOPS data. Thirdly, the spatial variance of GTLR-BiSAR data were thoroughly analyzed based on the range-Doppler (RD) geolocation functions. On the basis of a former analysis, we put forward the azimuth variance correction strategy and modiﬁed the conventional chirp scaling function to solve the range variance problem. Finally, we completed the derivation of the two-dimensional spectrum after the range chirp scaling. On the basis of spectrum expressions, we compensated for the quadratic and residue phase, and the azimuth compression was completed by SPECAN operation. In addition, we provide a ﬂow diagram to visually exhibit the processing procedures. At the end of this paper, the simulation and real data experiment results are presented to validate the effectiveness of the proposed algorithm.


Introduction
Bistatic synthetic aperture radar (BiSAR) has received increased interest throughout the world in the past decades due to its broad application prospects. Many researchers from different countries have carried out experiments with different platform pairs [1][2][3][4][5][6]. There appears to be a novel kind of bistatic remote sensing configuration whose transmitter is placed on geostationary satellite (GEO) with receivers placed on the low earth orbit satellite(LEO) [7][8][9][10], named the GTLR-BiSAR system. Compared with the monostatic case, the separation characteristic of GTLR-BiSAR has benefits, such as lower military vulnerability, saved costs, anti-interference, more abundant scattering information, and real-time monitoring, etc. [11,12]. These advantages make the GTLR-BiSAR system an appropriate spaceborne system for Earth observation missions with high-resolution, wideswath, high flexibility, and shorter revisit cycle characteristics [13,14]. However, the timevarying relative positions of the GEO transmitter and LEO receivers bring challenges to their methods were both proposed for the stripmap working mode, the azimuth Doppler aliasing problem in SS and TOPS modes will make their methods invalid. The geometry and TOPS working mode in [25] are more relevant to our research, where Sentinel-1 serves as a non-cooperative by using the TOPSSAR technique and the software-defined radio (SDR) hardware works as a fixed receiver. Facing the azimuth Doppler aliasing problem, Wang and Lu [26] adopted the azimuth multichannel to recover the unambiguous signal. Their method was not suitable for the single-channel receiver operating in SS or TOPS modes. In brief, the existing methods were applied on condition that the signal was not aliasing in azimuth time and the frequency domain. Hence, in this paper, we put forward the azimuth variance correction strategy and modified the conventional chirp scaling function to solve the range variance problem in the range-Doppler (RD) domain based on their research. Then the azimuth compression was completed in the azimuth frequency domain with the idea of spectral analysis (SPECAN) processing.
The remainder of this paper is organized as follows. Section 2 shows the imaging geometry and signal properties of the GTLR-BiSAR system; this is the basis of our proposed algorithm. Section 3 introduces the geometrical and mathematical methods to calculate the azimuth frequency modulation rate. The reference function is given to solve the azimuth Doppler aliasing problem. Section 4 describes our strategies to complete the azimuth and range variance corrections in the RD domain. Section 5 explicates the two-dimensional spectrum expression of the signal after the chirp scaling operation, and the azimuth compression is completed in the frequency domain. Section 6 shows the experiments of the simulated data. The correctness of our proposed processing algorithm has been verified. Furthermore, Section 7 introduces an equivalent real data experiment for demonstrating the proposed algorithm. Section 8 concludes the paper briefly.

GTLR-BiSAR Imaging Geometry and Signal Properties
In this section, the geometry of the GTLR-BiSAR system is described first. Then, the signal model is established and the signal expressions in different domains are deduced and analyzed in detail.

GTLR-BiSAR Imaging Geometry of Sliding Spotlight and TOPS Modes
As shown in Figure 1, the coordinate frame O-XYZ is utilized to explicate the geometry of the GTLR-BiSAR system. We assume the X-axis is the range direction of the illuminated scene. The Y-axis is the nadir of the LEO satellite tract, which is parallel to the flight direction. We assume it as the azimuth direction. The Z-axis is determined by the righthand Cartesian coordinate system, which is vertical to the illuminated scene plane and points to the Earth center inversely. It can be used to express the orbit height of the LEO satellite. The coordinate origin O is the starting nadir point of the LEO satellite.
The GEO transmitter is geostationary with the capability of monitoring the earth continuously. The main beam of the LEO antenna is steered to point to a fixed point called virtual rotation center, Figure 1a describes the SS mode of the LEO receiver, which is able to achieve high spatial resolution by extending the synthetic aperture time. Figure 1b describes the TOPS mode of the LEO receiver, which is able to provide excellent coverage performance by increasing the beam forward and backward steering angle. The black dashed lines L 1 and L 2 are parallel to the Y-axis, they are the far and near edge lines of the illuminated scene, respectively. The red dashed line L 3 is the center line of the swath. The red cross stars T 1 , T 2 and T 3 scatter along L 3 (T 3 is half red and half green because it is the intersection of L 3 and L 4 ). These targets along L 3 are characterized with the same Doppler FM rate because of the minimal range from these targets to the moving station, i.e., the LEO receiver, are all the same [23]. However, the bistatic ranges of these targets are different, their echoes are received in different range gates. The green cross stars T 3 , T 4 and T 5 scatter along L 4 . These targets are characterized with the same bistatic range (at a certain time). However, the FM rates of these targets confined in the same range gate differ from each other because their minimal ranges to the LEO receiver are different [23]. That is the main reason why the traditional algorithm fails to obtain a well-focused image, which is generally called the spatial variance problem in the BiSAR system. In Section 4, we will describe how to correct the spatial variance problem in detail.

GTLR-BiSAR Signal Model and Properties
In this subsection, we will deduce the signal expressions in different domains, including time domain, two-dimensional (2D) frequency domain (FD), and range-Doppler (RD) domain. Assume that a certain target point P is located in the illuminated scene, the minimal range from the target P to the LEO receiver is R r 0 p , and the range from the target P to the GEO transmitter is R T p . The effective radar velocity (see Appendix A for details) of the target P is V p , and the azimuth time is t a . It should be noted that all the velocity quantities in the rest of this paper refer to equivalent velocity. Then, the bistatic time-varying range from the target P to the GEO transmitter and LEO receiver is expressed as follows: Furthermore, the signal model of P can be written as: where K r is the FM rate of the chirp signal. f c is the carrier frequency. c is the velocity of light. τ is defined as the range time. W( t a T a ) and W( τ T τ ) are the azimuth and range time window. T a is the azimuth synthetic aperture time and T τ is the pulse duration. Ignoring unimportant amplitude window, we pay more attention to the phase terms in the following derivation.
Firstly, we perform the Fourier transform of Equation (2) to obtain the range spectrum by adopting the principle of the stationary phase (POSP) [27]: where f τ is the range frequency. Assuming that the integral phase of Equation (3) is θ, then By calculating the partial derivative of θ to τ and finding the zero point, we obtain the relationship between τ and f τ .
Let us bring τ (Equation (6)) into θ (Equation (4)), the range spectrum expression is obtained as follows: Secondly, we perform the Fourier transform of Equation (7) in the azimuth time to obtain the 2D frequency domain expression: where f a is the azimuth Doppler frequency. The integral phase θ of Equation (8) becomes: Similar to the former range operation, we obtain the relationship between azimuth time t a and Doppler frequency f a .
We set Z( f a , V p ) = , and rewrite Equation (14) in ascending power of f τ , successively: We obtain the concise 2D frequency spectrum expression. Finally, we perform inverse Fourier transform of Equation (16) in the range dimension to obtain the RD domain expression: where f s is the sample rate of the receiver. The integral phase θ of Equation (17) becomes: Similar to the former operations, we obtain the relationship between f τ and τ.
Submit f τ into Equation (18), the RD spectrum is given as follows: So far, the signal expressions in time, 2D frequency, and RD domains have been derived in detail. They are crucial for the design of the GTLR-BiSAR imaging algorithm for the SS and TOPS modes. We will introduce the specific processing procedures according to the signal properties in the following sections.

Azimuth Deramping Pre-Processing for GTLR-BiSAR
During the data acquisition of GTLR-BiSAR operating in SS or TOPS mode, the azimuth Doppler bandwidth of the echo data increase with the steering of the LEO receiver antenna's main beam. Under this circumstance, the azimuth Doppler bandwidth is much larger than the pulse repetition frequency (PRF), which results in the azimuth signal aliasing in the Doppler domain. In other words, the supporting area of the signal azimuth frequency spectrum for GTLR-BiSAR SS or TOPS modes is folded because of spatial undersampling. In order to process the echo data of GTLR-BiSAR SS or TOPS modes efficiently in the frequency domain, we should remove the azimuth aliasing first.
As shown in Figure 2, the coordinate frame O-XYZ is consistent with the previous definition in Figure 1. The black dashed lines L 1 and L 2 indicate the edges of the illuminated scene. The black dotted line L c stands for the center line of the illuminated scene. These lines are parallel to the flight direction of LEO receiver. The LEO receiver moves with an effective radar velocity V e f f , and the closest range from the virtual rotation center to the LEO receiver's trajectory is R re f . The range between L c and the LEO receiver's trajectory is R r0 . According to Figure 2, the equivalent projecting velocity of the scene center line is V c = V e f f /α geometrically. For the GTLR-BiSAR SS mode (Figure 2a), α = R re f /(R re f − R r0 ), and for the GTLR-BiSAR TOPS mode (Figure 2b), α = R re f /(R re f + R r0 ). We define α as the sliding factor, which is a critical parameter to determine the azimuth resolution and scene width for the SS and TOPS modes. To remove the azimuth aliasing, Carrara et al. [19], Lanari et al. [20], and Sun [21] have proposed the azimuth pre-processing methods for the SS or TOPS mode in monostatic SAR system. In this section, we further propose two kinds of uniform pre-processing methods, mathematically and geometrically, for the SS and TOPS modes in the GTLR-BiSAR system. The two methods are essentially the same. They are both utilized to calculate the azimuth frequency modulation rate K rot caused by the rotation of the LEO receiver main beam.
(1) The geometrical method to obtain K rot Using Figure 2, we can calculate the Doppler history f d_rot of the virtual rotation center according to the LEO receiver orbit track data, i.e., the position and velocity data of the LEO satellite. Assume that the azimuth duration of the LEO receiver orbit track data are t track , K rot is calculated as follows: where Max(·) and Min(·) indicate the maximal Doppler value and the minimal Doppler value, respectively.
(2) The mathematical method to obtain K rot According to the conventional monostatic SAR deramping operation in SS or TOPS modes, K rot is usually calculated by using the parameters R re f , λ and V e f f , as follows: For the GTLR-BiSAR system operating in SS or TOPS modes, it is worth noting that only the moving station, i.e., the LEO receiver, provides the azimuth Doppler contribution, so the constant coefficient 2 is absent in Equation (23) compared with the monostatic case.
Once we obtain K rot by either of the above mentioned methods, the reference function can be expressed as: Furthermore, the azimuth deramping pre-processing to remove the azimuth Doppler aliasing is achieved through a convolution operation with the above reference function: where ⊗ Azi denotes the azimuth convolution, and t a is the new azimuth time corresponding to the new azimuth frequency f a . It should be noted that ξ is a switch factor with the values of ±1. When ξ = 1/ − 1, the convolution operation is utilized to remove the Doppler aliasing of SS/TOPS mode. The convolution operation can be divided into two parts: the integration part and the phase compensation part. After the removal of azimuth aliasing, the azimuth time and frequency are updated so that we can process the echo data in the azimuth frequency domain efficiently. The detailed explanations of these new parameters and the convolution operation steps can be referred in [21], so we will not go into much detail in our paper. It should be noted that we still utilize t a and f a to express the updated azimuth time and frequency in order to be consistent with the previous and subsequent derivations in Sections 2, 4 and 5. We focus on the correction of spatial variance for GTLR-BiSAR system in the next section.

Spatial Variance Analysis and Correction for GTLR-BiSAR
In this section, we focus our attention on the spatial variance characteristics of GTLR-BiSAR system. The bistatic range-Doppler geolocation functions are introduced first. On this basis, the azimuth variance characteristic in the range-Doppler domain is expounded in detail. Additionally, we propose a novel method to solve the azimuth spatial variance problem effectively. Finally, the range variance characteristic in the RD domain is described in detail. Furthermore, we derive the chirp scaling (CS) function, especially for GTLR-BiSAR system to correct the range variance problem.

Range-Doppler Geolocation Functions for GTLR-BiSAR System
As shown in Figure 3, the bistatic signal propagation path starts from the GEO transmitter, it propagates from GEO transmitter to the target (the red arrow solid line points to the red star). Through back-scattering, the signal is received by the LEO receiver (the blue arrow solid line points to LEO satellite). The coordinate frame O-XYZ, dashed lines L 1 L 2 and the virtual rotation center are consistent with the previous definitions in Figure 1. We use S G and V G to represent the position and velocity vector of GEO transmitter. Correspondingly, S L and V L stand for the position and velocity vector of LEO receiver.
The position vector of the star target is assumed as S T . The descriptions of bistatic RD geolocation functions are not dependent on the operation mode, so we only give the geometry of GTLR-BiSAR system operating in the SS mode in Figure 3. According to Figure 3 above, we obtain S T by simultaneously solving the three geolocation equations, including the ellipsoidal equation, the Doppler equation, and the round-trip slant range equation [28,29]: where S T = S Tx , S Ty , S Tz , R G and R L stand for the distances from GEO transmitter and LEO receiver to the target, respectively. λ is the wavelength of transmitted signal. f dc is the Doppler center of the GTLR-BiSAR echo data, which is mainly determined by LEO receiver since GEO transmitter is geostationary. R e is the equatorial radius of Earth, and R q is the polar radius of Earth.It is obvious that the GTLR-BiSAR RD geolocation equations are nonlinear, Newton iterative method is adopted to obtain the target position (See Appendix B for details). According to Equation (26) and the Newton iterative method expatiated in Appendix B, we can obtain the target geolocations in the illuminated scene. We calculate the latitude and longitude coordinates of 25 point targets and mark them in Google Earth map with pushpins shown as Figure 4. The pushpins with the same color indicate the targets with the same minimal range Min(R L ) to the moving LEO receiver. Referring to the previous Equation (21), which shows the RD spectrum expression of a certain target P, we know that R T p and R r0p (Min(R L )) varying with the geolocation distribution of targets in the illuminated scene result in the spatial variance of the RD spectrum. The Doppler FM rate mainly depends on R r0p and the RCMs are relevant with both R T p and R r0p .
For the echo of GTLR-BiSAR in the RD domain, the targets confined in the same bistatic range gate differ from each other in R T p and R r0p . The spatial variance of R T p and R r0p causes the variance of Doppler FM rate and RCMs of the targets along the azimuth direction. As a result, the traditional RD and CS algorithms become invalid. Moreover, the existing azimuth perturbation function generated by nonlinear chirp scaling (NLCS) or local fit method [11,23] cannot be applied to equalize the targets' azimuth FM rates directly because the signal is aliasing in time azimuth domain after the deramping operation in Section 3. Therefore, we propose a novel method to solve the GTLR-BiSAR azimuth and range variance problem in Sections 4.2 and 4.3.

Azimuth Variance Analysis and Correction Strategy
As shown in Figure 4, the yellow pushpins along the azimuth direction represent the targets cluster 3, 8, 13, 18, 23. On one hand, the cluster is characterized with the same minimal range to the moving LEO receiver, denoted by R r0_re f . On the other hand, we suppose the distances from these targets to GEO transmitter as R T3 , R T8 , R T13 , R T18 and R T23 . It is obvious that these distances are different from each other. We suppose R T13 as the reference GEO distance R T_re f , if we can adjust the GEO distance differences and rearrange them with the same R T_re f . Then the signal azimuth variance will be removed and the traditional RD or CS method can be applied to process the azimuth variance corrected GTLR-BiSAR echo data. The azimuth variance correction schematic is shown in Figure 5, and we only give three target clusters for distinct illustration.   Figure 4, the supporting area of the targets cluster with the same R r0 spread to different bistatic range gates R bi . The vertical direction stands for the Doppler frequency f a , and f a_re f is the azimuth reference frequency. The discussion in this paper is assumed as the zero-Doppler, so f a_re f = 0. The horizontal direction represents the bistatic range gate, R bi_near , R bi_re f and R bi_ f ar are the near, reference and far bistatic range gate, R bi_re f = R T_re f + R r0_re f . The specific correction steps are shown in Figure 6.  Secondly, we multiply a linear phase function h i ( f r ) in the RD domain with the selected RD data (the red and blue dashed boxes in Figure 6). After inverse Fast Fourier transform (IFFT) operation in range direction, the RD data migrates to the desired positions (the right part of Figure 6). The migrations of targets 1-5 are R T i -R T (i+10) (i ∈ [1,5], i ∈ Z + ), recorded as ∆R T i . According to GEO satellite orbit tracks and target geolocations, the migrations are calculated and they change slightly along the range direction. Therefore, we choose ∆R T 3 to calculate ∆t 1 in Figure 6 to complete the signal migration. It should be noted that the phase compensation is needed after the signal migration. Since the phase is sensitive, each ∆R T i is utilized to calculate the compensation phase, and the compensation function is shown as follows: After the migration and compensation operations, the azimuth variances of targets 1-5 are removed. For the other targets, the similar operations are applied to complete the azimuth variance removals.
Through the azimuth variance corrections, the signal in RD domain becomes the the right part of Figure 5, which means that the target clusters with the same R r0 are characterized with the same R T . Referring to the previous Equation (21), we make a conclusion that the conventional monostatic processing strategy of RD and CS algorithms can be modified to handle the range variance problem. In order to avoid the blocking processing of RD algorithm, we prefer to modify CS algorithm to handle the range variance problem in next subsection.

Range Variance Analysis and Correction Strategy
After the correction of azimuth variance of RD data, the supporting area of the targets cluster with the same R r0 migrates to the desired positions in RD domain. Nevertheless, the supporting area of the targets cluster is characterized with different curvatures along the range direction shown in Figure 8. Facing the range variance problem, we firstly analyze the bulk, total and diff RCMs in detail. After that, we deduce the chirp scaling function for GTLR-BiSAR system to equalize the curves on the basis of the former RCMs analysis. Finally, the correction strategy of the range variance is completed. The vertical and horizontal directions of Figure 8 are defined as the azimuth frequency and the bistatic range gate, in accordance with the definitions in Figure 5. The horizontal dashed line f a and f a_re f are the Doppler frequency and azimuth reference frequency, the vertical black dotted lines represent the near R bi_near , middle R bi_re f , and far R bi_ f ar bistatic ranges. The solid green, red, and blue curves are utilized to express the RCMs of targets located in the near, middle and far bistatic range gates, respectively. The three red dotted curves represent the RCMs of targets located in the middle bistatic range gate. By shifting the red dotted curve to the near and far range gates, we distinguish the range variance characteristics visually. We assume the distance from vertical black dotted lines to the responding solid curves as RCM total . The distance from vertical black dotted lines to the red dotted curves is supposed as RCM bulk . The difference between RCM total and RCM bulk is defined as RCM di f f . According to Equation (21) in Section 2, we obtain the different RCM expressions as Equation (28). Take P as an example (the blue point in Figure 8), P is located in the far bistatic range curve. P c is located in the middle bistatic range curve. Connect P c with P, the black dashed line P c -P intersects the red dotted curve at P'. If we can shift P to P to complete the range variance removal, the middle and far curves are characterized with the same RCM bulk , the unified RCM correction can be performed to improve the processing efficiency.
The RCM di f f shift can be achieved by multiplying a chirp scaling function (CSF) [23] in monostatic SAR processing. On this basis, we deduce the CSF for the GTLR-BiSAR system. According to Equation (21) in Section 2, the target P locates in the range gate (in unit of time) as follows: set P c (the red point in Figure 8) as the original time, then we get: Express RCM di f f with τ as follows: Since V r changes slightly along the range direction, the third term of ∆τ (τ , f a ) can be ignored. Furthermore, we obtain the CSF for GTLR-BiSAR system by the following integral: Multiply the azimuth variance corrected signal with S CS (τ , f a ), the RD signal is characterized with the same RCM bulk , the unified RCM correction can be applied to remove the range variance problem for GTLR-BiSAR system.

Modified Chirp Scaling and Azimuth Compression for GTLR-BiSAR system
In this section, we firstly concentrate our attention on the 2D spectrum deduction of the signal after range chirp scaling operation in RD domain. Then, we complete the azimuth compression and focus the image in the azimuth frequency domain because the signal is aliasing in azimuth time domain after the deramping pre-processing. Finally, we provide a flow diagram to illustrate the entire process of our proposed algorithm visually.

Two-Dimensional Spectrum Expression after Chirp Scaling Operation
According to the range variance analysis in Section 4.3, we should multiply the RD signal (referring to Equation (21) in Section 2) with the chirp scaling function (referring to Equations (30) and (32) in Section 4.3) to remove the range variance: We adopt POSP again to deduce the 2D spectrum expression, and the integral phase θ(τ, f a , f τ ; P) is shown as follows: Calculate the partial derivative of θ to τ and find the zero point, we obtain the relationship between τ and f τ : Submit τ into Equation (34), we obtain the 2D frequency spectrum in ascending power of c successively: After proper derivations, Equation (37) is simply expressed as follows: The six exponential terms in the above formula are explained as follows: • The first term is related to azimuth Doppler frequency modulation. It is a function of GEO distance R T p and minimal LEO distance R r0p , which will be compensated during azimuth processing. We assume this term as H az ; • The second term is a linear phase in range frequency domain, it represents the peak position after range compression; • The third term stands for the bulk RCM after chirp scaling operation. The RCM correction of the whole scene is completed by compensating this phase, recorded as H RCM ; • The fourth term means the range modulation after chirp scaling operation. It is a quadratic function of f τ . The range compression is achieved by compensating this phase, assumed as H ran ; • The last two terms is an residual phase caused by chirp scaling operation, they are functions of bistatic range R bi and azimuth frequency f a . We define this term as H RP . They should be removed during azimuth processing.
In the traditional monostatic CS algorithm processing for stripmap workmode, there are also similar exponential phase terms to be compensated. Different from GTLR-BiSAR system operating in SS or TOPS modes, there is no azimuth aliasing problem for the former system. Therefore, we should focus the image in the azimuth frequency domain after we complete the phase compensation operations in 2D frequency domain [21].

Azimuth Compression in Azimuth Frequency Domain
For the traditional monostatic stripmap workmode, the azimuth compression is achieved by compensating the azimuth nonlinear phase in RD domain and changing the data through azimuth IFFT. Then the azimuth compression is completed in time domain. However, the data of GTLR-BiSAR operating in sliding spotlight or TOPS mode is aliasing in azimuth time domain after the deramping pre-processing ( Section 2). Combined with the idea of SPECAN processing, we focus the point target in the azimuth frequency domain as follows.
After the range compression, bulk RCM correction, residual phase compensation and azimuth phase removal ( referring to Equation (38) in the former subsection), we construct a transfer function as follows: Then we obtain the signal expressions: where * represents the conjugate operation. The S t f ( f a ; P) ( f τ is removed in the operation above) is no longer aliased in azimuth time domain, so we convert it to azimuth time domain by azimuth IFFT, and the signal becomes: Compensate the above phase in time domain and transfer it to azimuth frequency domain, we finally obtain the focused image.

The Entire Flow Diagram of The Proposed Algorithm
In order to promote better understanding of the entire processing procedures visually, we provide a flow diagram (Figure 9) where the azimuth deramping pre-processing, spatial variance correction, and azimuth compression are shaded pink in individual blocks.

Simulation Experiments and Results
In this section, we will introduce the simulation conditions and experiment results to validate our proposed algorithm. Taking the GTLR-BiSAR sliding spotlight mode as an example, we will elaborate on the entire simulation experiment and processing details, shown as Figure 10. Table 1 describes the experiment parameters of GEO and LEO satellites. The satellite orbit data are generated by satellite tool kit (STK) software and the orbit data are interpolated according to pulse repetition interval (PRI). The parameters of GEO transmitter and LEO receiver are also listed, which are designed according to the desired resolution. The illuminated scene size is determined by the antenna main lobes of GEO and LEO antennas. The point targets are geolocated as shown in Figure 4. Based on the simulation parameters above, we give the results of simulation and processing and make appropriate analysis.  The typical imaging performance indexes of the targets focused by the proposed method are shown in Figure 11 and Table 2. It is worth noting that the range resolution is related to the bandwidth of the transmitted signal. In order to reduce the amount of simulation data, the resolution is designed as 0.8 m (azimuth) × 1.5 m (range) theoretically. The point targets of GTLR-BiSAR system operating in sliding spotlight mode is shown in Figure 4. The 5 × 5 point targets are scattered in the illuminated scene area with the size of 5 km (azimuth) × 5 km (range).  Table 2. PSLR, ISLR, and IRW of selected targets of SS mode in Figure 11 by the proposed method.

Target
Range Azimuth In addition, the convention NLCS method is utilized to complete the imaging operation. We take the scene center point , i.e., the point target 13 as an example, then we compare the imaging results of the convention NLCS method and the proposed method. The comparison results show that the NLCS method cannot be applied in the GTLR-BiSAR sliding spotlight case, because the signal is azimuth aliasing in time domain after the azimuth deramping pre-processing in Section 3, while the NLCS method is suitable for the case that the signal is not aliasing in both time and frequency domain. As a result, the point target is defocused seriously and the image is blurring, shown in Figure 12.

PSLR (dB) ISLR (dB) IRW (m) PSLR (dB) ISLR (dB) IRW (m)
After the experiment of GTLR-BiSAR sliding spotlight mode, we design the parameters and carry out simulation experiments and data processing of GTLR-BiSAR TOPS mode. The simulation parameters and data processing results are are given directly as follows. Since the processing procedure is similar to the sliding spot mode, the focused and evaluation results are given directly. Similarly, the NLCS method cannot be used in this case. Table 3 describes the experiment parameters of GTLR-BiSAR TOPS mode. Compared with Table 1 the sliding factor is different because of the geometry difference. The foused results are shown in Figure 13 and Table 4. The NLCS result is shown in Figure 14. Moreover, the transmitted bandwidth is much smaller because GTLR-BiSAR TOPS mode is characterized with large-scale imaging and lower resolution. The resolution is designed as 10 m (azimuth) × 10 m (range) theoretically. The 5 × 5 point targets are scattered in the illuminated scene area with the size of 80 km (azimuth) × 40 km (range).   Figure 13. The imaging results of point targets for GTLR-BiSAR TOPS mode. Table 4. PSLR, ISLR, and IRW of selected targets of TOPS mode in Figure 13 by the proposed method.  In the practical application of SAR image formulation, there exists parameters measurement errors, such as the velocity error and slant range error. In order to reduce the image defuse, we can integrate some autofocus algorithms [30] with our proposed approach to improve the focusing accuracy. Moreover, we apply the traditional backprojection (BP) method to focus the center target. The result emphasizes the significance of our proposed approach. The PSLR, ISLR and IRW of selected target by time domain method are basically consistent with the frequency method we put forward. However, it takes 3.6798 s to focus a single target, which is an unbearable time consumption for scene imaging. So, we focus our attention on the frequency domain method and apply it in the scene imaging of wide swath.

Equivalent Real Data BiSAR Experiment Based On BeiDou Navigation Satellite
In order to further verify our frequency domain imaging algorithm, we designed and built a satellite-ground SAR system to carry out equivalence experiment. The real data experiment is organized as Figure 15. In the experiment, BeiDou navigation satellite with inclined geosynchronous orbit (IGSO) serves as the transmitter operating in the L-band. The transmission signal adopts binary phase shift keying (BPSK) modulation mode with the bandwidth of 20.46 MHz. In our experiment, we found that the reflected signal was too weak to obtain a SAR image. Therefore, we set a repeater in the illuminated scene to amplify the reflected signal to improve the signal-to-noise ratio (SNR). As shown in Figure 15a, the coordinate system O − XYZ is consistent with the definition in Section 2.1.
In the scene flat, the antenna R receives the enhanced signal from the repeater. A few tens of meters away from antenna R, there exists antenna D which receives the direct signal from BeiDou navigation satellite for synchronization. The received echo is transferred to the computer for focusing processing.
It should be noted that our GTLR-BiSAR focusing algorithm is proposed for the bistatic configuration of geostationary transmitter with LEO receivers, where the emitter is static and the LEO receiver is moving. Our developed formalism can be applied to such a geometry where BeiDou navigation IGSO satellite is a moving emitter source with a static ground receiver. Therefore, our proposed method can be demonstrated by the real data experimental results to a certain extent. Figure 15b shows our experiment conducted on the roof. Some main parameters are given in the following Table 5.  Different from the traditional spaceborne SAR satellites, BeiDou IGSO satellite transmits BPSK continuous navigation signal. Firstly, we should transform the continuous signal to segmentation form according to the coarse/acquisition (C/A) code. Then we obtain the two dimensional echo matrix shown in Figure 16a. Then we apply our proposed method to focus the image, the result in Figure 16b-d validate the effectiveness of our algorithm. The theoretical range resolution is 19.40 m, and the result of our method is 20.99 m. The theoretical azimuth resolution is 4.55 m, and the result of our method is 4.62 m. The non-ideal azimuth pulse response is caused by the system errors. It should be noted that the range compression is operated according to the BPSK code, which is different from the chirp signal. Essentially, they are convolution operations to realize the signal autocorrelation.

Conclusions
In a GEO-LEO BiSAR system with the characteristics of complex geometry and obvious orbit difference, the signal of the illumined area is spatial variant because the relative position of the transmitter and the receiver changes with time. By conventional NLCS method, azimuth perturbation function is introduced by local fit method on condition that the signal is not aliasing in time and frequency domain. However, the method becomes invalid for the bistatic SAR configuration of geostationary transmitter with LEO receivers working in sliding spotlight or TOPS modes.
In this paper, a unified algorithm for the GTLR BiSAR sliding spotlight and TOPS modes data processing is proposed. Firstly, the GTLR-BiSAR imaging geometry and signal properties are introduced. Then, the azimuth Doppler aliasing problem is solved by deramping operation. Furthermore, the spatial variance problem is corrected by our proposed strategies. In addition, the modified chirp scaling function is proposed to complete the RCM correction and the azimuth compression is operated in frequency domain. Finally, simulation and real data experiments are carried out to validate our algorithm. In the future, interferometry processing methods will be studied in the GTLR-BiSAR system for the generation of digital elevation model. Acknowledgments: The authors would like to thank anonymous reviewers for their helpful comments and suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In this appendix, we will introduce the least squares method to calculate the equivalent velocity of the LEO satellite. The fitted polynomial of the range from the satellite to the target is supposed as y . a 0 and a 1 are the fitting coefficients of the constant and first order term, respectively. a k is the fitting coefficient of the k th (k nd , k rd , k ∈ N, N is the natural number) order term. Then, y is expressed as follows: Assume that the curve to be fitted is y, then the sum of distances between y and y can be expressed as follows: [y i − (a 0 + a 1 x i + ... + a k (x i ) k )] 2 (A2) Find the minimum value of the above variance, and the partial derivatives of a 0 , ...,a k are sequentially calculated as 0: After derivations of the above formula, we obtain the following equations: Express the above equations in the form of matrix: We transform the first item of the matrix equation to the Vandermonde form: x 1 x 2 · · · x n . . . . . . . . . . . . , then the above matrix equation can be simplified as X XA = X Y. The coefficient vector A is obtained by A = (X X) −1 X Y, and the effective radar velocity is calculated according to the coefficients.

Appendix B
In this appendix, we will introduce the Newton iterative method to solve the nonlinear GTLR-BiSAR range-Doppler geolocation equations. Assuming that the position and velocity vectors of LEO receiver are S L = (S Lx , S Ly , S Lz ) and V L = (V Lx , V Ly , V Lz ), separately. Correspondingly, S G = (S Gx , S Gy , S Gz ) and V G = (V Gx , V Gy , V Gz ) stand for the position and velocity vectors of the GEO transmitter, respectively. Then we suppose the range vector from the LEO receiver to the target as S LT = S T − S L = (x, y, z), where x, y and z are the three-axis components of the vector. Moreover, the range vector from the GEO transmitter to the target is assumed as S GT = S T − S G = S T − S L + S L − S G = (x + S Lx − S Gx , y + S Ly − S Gy , z + S Lz − S Gz ) = (x + ∆x, y + ∆y, z + ∆z), where we introduce S LT and the vector from the GEO transmitter to the LEO receiver. ∆x, ∆y and ∆z are the three-axis components of the vector from the GEO transmitter to the LEO receiver. Finally, we can change the GTLR-BiSAR range-Doppler geolocation equations into the following functions: R q 2 f 2 (x, y, z) = xV Lx + yV Ly + zV Lz (x + ∆x) 2 + (y + ∆y) 2 + (z + ∆z) 2 + (x + ∆x)V Gx + (y + ∆y)V Gy + (z + ∆z)V Gz x 2 + y 2 + z 2 −λ f dc x 2 + y 2 + z 2 (x + ∆x) 2 + (y + ∆y) 2 + (z + ∆z) 2 f 3 (x, y, z) = R R + R L − x 2 + y 2 + z 2 − (x + ∆x) 2 + (y + ∆y) 2 + (z + ∆z) 2 (A7) Supposing that X =   x y z   , f (X) =   f 1 (x, y, z) f 2 (x, y, z) f 3 (x, y, z)   and X k is an approximate zero point of f (X). f (X) is approximated by the linear part of the Taylor expansions (the first two terms) at this zero point: