Improved GNSS-Based Bistatic SAR Using Multi-Satellites Fusion: Analysis and Experimental Demonstration

The Global Navigation Satellite System (GNSS)-based Bistatic Synthetic Aperture Radar (SAR) is getting more and more attention in remote sensing for its all-weather and real-time global observation capability. Its low range resolution results from the narrow signal bandwidth limits in its development. The configuration difference caused by the illumination angle and movement direction of the different satellites makes it possible to improve resolution by multi-satellite fusion. However, this also introduces new problems with the resolution-enhancing efficiency and increased computation brought about by the fusion. In this paper, we aim at effectively improving the resolution of the multi-satellite fusion system. To this purpose, firstly, the Point Spread Function (PSF) of the multi-satellite fusion system is analyzed, and focusing on the relationship between the fusion resolution and the geometric configuration and the number of satellites. Numerical simulation results show that, compared with multi-satellite fusion, dual-satellite fusion is a combination with higher resolution enhancement efficiency. Secondly, a method for dual-satellite fusion imaging based on optimized satellite selection is proposed. With the greedy algorithm, the selection is divided into two steps: in the first step, according to geometry configuration, the single-satellite with the optimal 2-D resolution is selected as the reference satellite; in the second step, the angles between the azimuthal vector of the reference satellite and the azimuthal vector of the other satellites were calculated by the traversal method, the satellite corresponding to the intersection angle which is closest to 90° is selected as the auxiliary satellite. The fused image was obtained by non-coherent addition of the images generated by the reference satellite and the auxiliary satellite, respectively. Finally, the GPS L1 real orbit multi-target simulation and experimental validation were conducted, respectively. The simulation results show that the 2-D resolution of the images produced by our proposed method is globally optimal 15 times and suboptimal 8 times out of 24 data sets. The experimental results show that the 2-D resolution of our proposed method is optimal in the scene, and the area of the resolution unit is reduced by 70.1% compared to the single-satellite’s images. In the experiment, there are three navigation satellites for imaging, the time taken to the proposed method was 66.6% that of the traversal method. Simulations and experiments fully demonstrate the feasibility of the method.


Introduction
Apart from providing the global positioning and timing services, the Global Navigation Satellite System GNSS could be used as the transmitters in bistatic radar to explore the earth parameters and detect targets all-day and all-weather. This technology is known as GNSS-Reflectometry (GNSS-R) which has the advantage of low-cost, global coverage because of no transmitters and massive GNSS satellites. At present, GNSS-R has been used to measure sea-surface level [1,2], retrieve wind speed [3,4], estimate soil moisture and biomass [5][6][7][8], detect oil spill [9] et al. Additionally, the imaging capacity of the GNSS-R signal has been also taken into account with the development of technologies.
GNSS-based Bistatic SAR is one emerging application that combines GNSS-R and bistatic SAR imaging technology. The combination of the two technologies can be selected at the signal level, navigation satellites are used as the non-cooperative illuminator, and the bistatic SAR mode is used for echo data collection and processing to generate imaging. The feasibility of the GNSS-based Bistatic SAR has been demonstrated with GPS [10,11], GLONASS [12], Beidou [13], and Galileo [14] as illuminators. Furthermore, a series of problems encountered during the research process were studied in depth, such as power budget [15], synchronization [16], and imaging algorithms [17]. This technology has been demonstrated to be meaningful for detecting [18], recognizing targets [19], and surface change monitoring [20,21]. Moreover, the number of GNSS satellites illuminating the same area on the ground from different aspect angles is about 20-30 for GNSS constellations. This factor makes it possible to select the best bistatic geometry or combine multiple satellites for observation. For the radar itself, it can improve resolution performance and reduce shadow area. In terms of remote sensing applications, it can effectively improve the probability of detection [22] and the accuracy of surface deformation monitoring [23].
GNSS signals are not designed for radar or remote sensing, prompting the introduction of several problems in radar and remote sensing applications. First of all, the power budget issue is affected by the transmitting power of navigation satellites and the orbit height (navigation satellites are mainly medium earth satellites: MEO, Medium Earth Orbit), resulting in very low power density near the surface of the earth. Secondly, the navigation signal bandwidth is relatively narrow to the traditional synthetic aperture radar signal, the range resolution is limited to the signal bandwidth. For example, the bandwidth of the GPS L1 signal is 2.046 MHz, which provides a 3 dB range resolution about 87.9 m in an equivalent monostatic model. With the further development of the GNSS systems, it provides a wider signal bandwidth. The GPS L5, Beidou B3, GALILEO E5a and E5b signals provide a 3 dB range resolution about 8.79 m. Since GNSS-based SAR is in the category of bistatic, the geometric configuration may greatly deteriorate the spatial resolution of the system [24], Therefore, the application field of this technology is limited.
To tackle the aforementioned problem, researchers have proposed different methods to improve imaging resolution. A sufficiently high signal-to-noise ratio (SNR) can be obtained by increasing the synthetic aperture time, which can be a few minutes or even longer for a fixed receiver [25]. Due to the orbit curvature of the navigation satellite, the long dwell time can achieve a substantial increase in the spatial resolution of the system [26]. One or more receivers receive multiple navigation satellites signals to generate different SAR images of the same target area, and the images can be fused in a coherently or incoherently way to improve the spatial resolution, it has been verified the feasibility and performance by experiments [27][28][29]. A large number of satellites with different views can be used to enhance the information space of the images, the feature value can be extracted from each image, and then the feature-level combination can be used to significantly enhance the target information [30], conducive to feature extraction, feature recognition, and classification of the target area. In [31], using the B3 signal of the BeiDou navigation system to study the segmentation method of the sensing area, the method is based on the changes of the radar cross section (RCS) of different targets under different bistatic geometric configurations, and is verified by experiments. With Galileo E5 echo imaging, the entire E5 channel can be stitched together using a spectral equalization technique to increase the 3 dB range resolution to 1.758 m, but this technique is limited to alternate binary offset carrier (AltBOC) signal and introduces essential reduction in signal-noise ratio (SNR) [14].
It has been proved that the resolution of the system can be improved by using space separated multi-satellite and multi-station imaging, however, these methods introduce new problems. The first is the increased computation, and the second is the resolution improvement efficiency problem caused by multi-satellite fusion. In this paper, to solve this problem, a multi-satellite fusion method with improved imaging resolution is proposed based on an optimized satellite selection strategy. We consider a ground-based stationary receiver configuration, but without loss of generality, the method proposed can be used in others topologies. The processing rules for signals from multiple satellites are combined based on multiple parameters (the transmitted signal, the configuration of the transceiver station, and the synthetic aperture time) defined by Multi-Satellite Fusion Point Spread Function (MPSF). In this paper, we focus on the influence of transceiver station configuration on the fusion resolution.
We established an accurate geometric configuration model and evaluated the influence of geometric configuration on the fusion resolution in the form of numerical simulation. Simulation results show that multi-satellite fusion can not only improve system resolution but also potentially reduce it, which is mainly affected by the geometric configuration of multiple satellites. The resolution improvement efficiency decreases with increasing number of satellites, and the geometric configuration requirements increase as the number of satellites increases. Additionally, the simulation reveals that dual-satellite fusion is a better cost-effective option because of the increased resolution efficiency. A fusion imaging method based on optimized satellite selection is proposed to realize the effective spatial resolution improvement of a multi-satellite system.
To demonstrate the feasibility of the proposed imaging algorithm, point target simulations were performed using real GNSS orbit data. The GPS L1 satellite is chosen as the illuminator. A total of 24 sets of data were collected in one day, and the proposed method was compared with the traditional traversal fusion method. The proposed method obtained optimal resolution images 15 times and suboptimal images 8 times, demonstrating the feasibility of the proposed algorithm.
An experiment is also conducted to further confirm the feasibility of the proposed algorithm. The GPS L1 satellite is chosen as the illuminator. The Beihang University Gymnasium was selected as the imaging target because of the large RCS. The collected echo data are processed using the traversal fusion method and the proposed algorithm, respectively. The proposed algorithm produces the optimal resolution image of the experimental scene, the time taken is 66.6% of the traversal fusion method, demonstrating the feasibility of the proposed algorithm.
The structure of this article is as follows: numerical analysis of PSF and MPSF are performed in Section 2, a fusion imaging method based on optimize satellite selection is introduced in Section 3, the simulation is carried out through the real navigation satellite orbit data, and the simulation results are analyzed in Section 4, the proposed method is verified and analyzed experimentally in Section 5, finally we summarize our conclusions in Section 6.

Generalized Ambiguity Function
In bistatic SAR, the ambiguity function (AF) of bistatic radar is not only related to the transmitted signal but also related to factors such as the bistatic topology structure. According to the research of M. Cherniakov and Tao Zeng et al. [32], the AF was extended to the GNSS-based Bistatic radar, thus the generalized ambiguity function (GAF) could be expressed as: where s A (t, u) and s B (t, u) are the GNSS baseband echo signals by a point target with unit radar cross section and located at the ground points A and B. O is center of the scene, → OA and → OB are the vectors from point O to point A and B, respectively. The variable u is called fast time, constant during the pseudo-random code period T (1 ms long for the GPS L1(C/A-code)), t denote the slow-time, set t = nT + u. K is used as normalization parameter such as χ ( → OA, → OA) = 1. In the article [33], the GAF simplified as follows: Sensors 2020, 20, 7119 is the matched filter output of the ranging signal, since navigation satellites emit continuous pseudo-random signals, the output of the matched filter is a triangular wave, which can be expressed as: where B is the navigation signal bandwidth, P 0 is power, tri(·) is triangular function, specifying that the range resolution. The delay τ d is related to the geometric configuration, which can be expressed as: where → β (u) is the bisector of the bistatic angle, c is the speed of light, u ta is the value of u evaluated at the midpoint of T (u ∈ [u ta − T/2, u ta + T/2]). m(u) represents the ratio of signal received power to transmit power, M( f ) is the inverse Fourier transform of m(u), M( f d ) is given by: sin c(·) specifies azimuth resolution. f d is Doppler frequency, which is generated in two parts: the motion of the receiver and the motion of the navigation satellite, which can be expressed as: where f is the navigation signal carrier frequency, λ is wavelength of navigation signal.

Resolution Analysis Parameter
We consider a ground-based stationary receiver collecting the signals emitted from a GNSS transmitter and reflected by a stationary scene. Figure 1 shows the geometric configuration, where the navigation satellite S is located on the OXYZ plane, the sub-satellite point S is located on the negative semi-axis of the Y-axis, R is a ground-based stationary receiver, and R is the ground mapping point, the satellite velocity → V T points to the positive direction of the X-axis, θ T x and θ R x is the ground incident angles of the navigation satellite and the receiver relative to the center O of the target area, the observation angle of the navigation satellite φ is the angle with where a r is a factor accounting the shape of g(t); because the signal transmitted by the navigation satellite is a pseudo-random code, the output of the matched filter is a triangular wave, for its 3dB attenuation, a r = 0.586. a a is a factor accounting the shape of M( f ); between the synthetic aperture time, the power of the received signal can be approximated as a rectangular function, so the output of the matched filter is a sin c function, which attenuates 3 dB, a a = 0.886. → β is the angular bisector vector of the satellite and receiver position, ω is the equivalent bistatic angular speed. The bistatic angle can be expressed as: Sensors 2020, 20, 7119

PSF Numerical Analytical
To more intuitively evaluate the two-dimensional (2-D) capability of GNSS-based Bistatic SAR, 2-D resolution unit was proposed in [32] to analyze the resolution capability of bistatic SAR, using the 3dB approximate elliptical area of its Point Spread Function (PSF) to characterize it [33]. The (PSF) is the result of (2) being mapped onto the ground plane, it is an approximately elliptical analytic unit. The PSF function is given by: The resolution ellipse is computed as: Formula (12) shows different bistatic PSF parameters are affected by the following factors: (1) synthetic aperture time T, (2) signal bandwidth B, (3) geometric configuration of bistatic. As shown in Figure 1, the influence of geometrical configuration on PSF is mainly considered, determined by a is the angle between the range and azimuth vector direction with the ground plane.
Since a ground-based stationary receiver, the imaging area will be relatively small, take the center O of the target area as a reference to calculate the ground-to-ground resolution and azimuth resolution. The unit vectors of range and azimuth projected to the ground are respectively:

PSF Numerical Analytical
To more intuitively evaluate the two-dimensional (2-D) capability of GNSS-based Bistatic SAR, 2-D resolution unit was proposed in [32] to analyze the resolution capability of bistatic SAR, using the 3dB approximate elliptical area of its Point Spread Function (PSF) to characterize it [33]. The (PSF) is the result of (2) being mapped onto the ground plane, it is an approximately elliptical analytic unit. The PSF function is given by: The resolution ellipse is computed as: Sensors 2020, 20, 7119 6 of 22 Formula (12) shows different bistatic PSF parameters are affected by the following factors: (1) synthetic aperture time T, (2) signal bandwidth B, (3) geometric configuration of bistatic. As shown in Figure 1, the influence of geometrical configuration on PSF is mainly considered, determined by navigation satellite position information (θ T x , θ R x and φ) and the velocity → V of the navigation satellite. Based on these 4 parameters we build accurate numerical simulation scenarios, assuming that the GPS L1 satellite is moving in a straight line at a constant speed. The specific parameters are shown in Table 1, and the simulation scene as shown in Figure 2, the results of PSFs can be evaluated using both numerical [33] and analytical methods, which can be expressed as A 1 (ϕ, γ, α), where α is the complementary angle of θ T x , ϕ is the supplementary angle of φ, θ R x ,γ is the angle with the satellite velocity → V s and the positive direction of the X-axis.  Table 1, and the simulation scene as shown in Figure 2, the results of PSFs can be evaluated using both numerical [33] and analytical methods, which can be expressed as 1 ϕ γ α    As shown in Figure 3a, the simulated optimal PSF A 1 (180 • , 90 • , 45 • ), η is the resolution ellipse's orientation, which is the function of the  Figure 4 shows the distribution of 2-D resolution units area produced to follow the change of ϕ and γ, the highlighted part shows that the area of the 2-D resolution unit is very large, indicating that the 2-D resolution capability extremely deteriorates under this geometric configuration.  Figure 4 shows the distribution of 2-D resolution units area produced to follow the change of ϕ and γ , the highlighted part shows that the area of the 2-D resolution unit is very large, indicating that the 2-D resolution capability extremely deteriorates under this geometric configuration.

MPSF with Non-Coherent Summation
The multi-satellite fusion system is a typical multi-static SAR system. The AF can be regarded as non-coherent summation of the AF of multiple bistatic systems [33]. The GAF of the multi-satellite fusion system can be expressed as:  Figure 4 shows the distribution of 2-D resolution units area produced to follow the change and γ , the highlighted part shows that the area of the 2-D resolution unit is very large, indicatin at the 2-D resolution capability extremely deteriorates under this geometric configuration.

. MPSF with Non-Coherent Summation
The multi-satellite fusion system is a typical multi-static SAR system. The AF can be regarde non-coherent summation of the AF of multiple bistatic systems [33]. The GAF of the multi-satelli sion system can be expressed as:

MPSF with Non-Coherent Summation
The multi-satellite fusion system is a typical multi-static SAR system. The AF can be regarded as non-coherent summation of the AF of multiple bistatic systems [33]. The GAF of the multi-satellite fusion system can be expressed as: m is the number of satellites participating in the fusion, B n is bandwidth of GNSS signal, P n is power, T n is synthetic aperture time, and → β n is the unit vector along the bisecting line of the bistatic angle of the system. It can be seen from the above Formula (14) that the GAF of the multi-satellite fusion system is related to the form of the transmitted signal, the configuration of the transceiver station, and the synthetic aperture time. In [28], the MPSF was defined as: The MPSF is a noncoherent combination of the individual PSFs, which resolution unit area is the intersection of the single-satellite PSFs. Using the 3dB projection area of the MPSF on the XOY plane as the 2-D resolution unit. Based on the radar parameters in Table 1 and the simulation scene as shown in Figure 2, we can be evaluated using both numerical and analytical methods. The results of MPSF can be approximated as Assume m = 2, the quantitative analysis of the MPSF are carried out, assume Figure 5a, its 2-D resolution area is 743 m 2 , compared with the A 1 (180 • , 90 • , 45 • ), the area is reduced by 47.2%. The 2-D resolution units of MPSF (m = 2) are shown in Figure 6a, Table 2 shows the relevant parameters of the 2-D resolution unit of the optimal MPSF. The value of the highlighted part in Figure 6a reveal that the area of the fused 2-D resolution unit has exceeded 1390 m 2 , indicating that the fusion has not improved the 2-D resolution, but worsened it, meaning that the fusion under this geometric configuration is meaningless. As shown in Figure 7b, compared to the optimal resolution (1390 m 2 ) of the single-satellite PSF, about 76.7% 2-D resolutions of MPSF (m = 2) had improved in varying degree, and about 23.3% 2-D resolutions were deteriorated.
Assume m = 3, the quantitative analysis of the MPSF are carried out, assume , the area of 2-D resolution units are shown in Figure 6b. Table 2 shows the relevant parameters of the 2-D resolution unit of the optimal MPSF (m = 3), the optimal MPSF is shown in Figure 5b, its 2-D resolution area is 616 m 2 , compared with the , the area is reduced by 55.7%. It shows that the fusion of three-satellite can obtain better resolution performance than the fusion of dual-satellite. Compared with the dual-fusion optimal solution (743 m 2 ), only about 30% 2-D resolution of MPSF (m = 3) has been improved in varying degrees, and about 70% 2-D resolutions were deteriorated. The results show that the threesatellite fusion has higher requirements for the geometric configuration of the satellites participating in the fusion than the dual-satellite fusion, and the complexity of satellite selection is further improved, but the performance improvement efficiency is reduced.

MPSF Numerical Simulation (m ≥ 3)
Assume m = 3, the quantitative analysis of the MPSF are carried out, assume A 1 (180 • , 90 • , 45 • ) + A 2 (100 • , 50 • , 45 • ) + A 3 (ϕ 3 , γ 3 , 45 • ), the area of 2-D resolution units are shown in Figure 6b. Table 2 shows the relevant parameters of the 2-D resolution unit of the optimal MPSF (m = 3), the optimal MPSF is shown in Figure 5b, its 2-D resolution area is 616 m 2 , compared with the A 1 (180 • , 90 • , 45 • ), the area is reduced by 55.7%. It shows that the fusion of three-satellite can obtain better resolution performance than the fusion of dual-satellite. Compared with the dual-fusion optimal solution (743 m 2 ), only about 30% 2-D resolution of MPSF (m = 3) has been improved in varying degrees, and about 70% 2-D resolutions were deteriorated. The results show that the three-satellite fusion has higher requirements for the geometric configuration of the satellites participating in the fusion than the dual-satellite fusion, and the complexity of satellite selection is further improved, but the performance improvement efficiency is reduced.
The MPSF results of the simulation were analyzed (assume m ≤ 9). As shown in Figure 5, the MPSF with a well-defined central response at the target position is surrounded by a series of side-lope. With the increase of m, the number of side-lopes increase, and the peak side-lobe ratio (PSLR) decreases. That will cause the image SNR to drop and require more post-processing work. As shown in Figure 7a, the optimal 2-D resolution unit area of the MPSF gradually decreases as m increases, indicating that the resolution improvement as the number of fusion satellites increases. However, the improved efficiency of resolution will be significantly reduced. The simulation results were shown in Figure 7b show that as m increases, compared with the optimal solution of m−1, the ratio of combinations with improved resolution decreases, while the ratio of combinations with decreased resolution increases. This means that as m increases, the geometric configuration requirements of the satellite combination with improved resolution will also increase, and the complexity of satellite selection will also increase. In this paper, we considered the following four aspects: resolution improvement efficiency, satellite selection complexity, introduce noise, and computation quantity. In summary that dual-satellite fusion is the cost-effective choice.

Improved Dual-Satellite Fusion Imaging Method
With the development of multiple navigation systems (GPS, Beidou-2, GLONASS, GALILEO), the average number of visible navigation satellites at the same place on the earth's surface will rise from 10 to about 40. Therefore, the observation information of GNSS-based Bistatic SAR systems will increase significantly and reduce the revisit period. However, the complexity of satellite selection and computation will increases follow the satellites number increase. According to the Section 2 of the simulation analysis, this section proposes an optimized dual-satellite fusion imaging method, obtaining the optimal 2-D resolution image of the dual-satellite fusion in the scene. To test the performance of this method, the method has been implemented in MATLAB R2018b (MathWorks, Natick, MA, USA) as an interpreted language, and because it is widely used in research and academia. The improved dual-satellite fusion imaging method is divided into two parts. The first part is the satellite selection strategy, the implementation of which is shown in the blue box in Figure 8 and will be explained in detail in Section 3.1; as shown in the red box in Figure 8, the second part is the fusion imaging algorithm, which we will elaborate on in detail in Section 3.2.
a. The improved dual-satellite fusion imaging method is divided into two parts. T e satellite selection strategy, the implementation of which is shown in the blue box in ill be explained in detail in Section 3.1; as shown in the red box in Figure 8, the second n imaging algorithm, which we will elaborate on in detail in Section 3.2.

Satellite Selection Strategy
For obtaining an optimal dual-satellite fusion resolution image, it is necessary to select the best combination from a large number of visible satellites for fusion imaging. To reduce computation and obtain better imaging resolution, a satellite selection method based on a greedy algorithm is proposed. The proposed method linearly approximates the best 2-D resolution of the dual-satellite fusion by successive greedy selection of two satellites in an iterative way. The specific steps are shown in Figure 8, Algorithm 1 gives the specific implementation of the pseudo code.

•
Step 1: Select a subset of navigation satellites for imaging Before imaging, using the navigation satellite ephemeris data to roughly determine the visible satellites distribution at the receiver's location, and build the geometric configuration model of GNSS based Bistatic SAR system based on the location information of the target area. Select a subset of navigation satellites for imaging based on the geometric model. In consideration of two important factors related to image quality, the first is the imaging resolution; the second is the signal-to-noise ratio (SNR) of the image, which is related to the received signal power. Therefore, the image quality is related to the geometric configuration of the GNSS-based bistatic SAR. The resolution can be determined by the Formula (7) in the range and azimuth direction. Combined with the power analysis of the navigation signal backscatter [34], the subset of navigation satellites for imaging can be constrained by constraining its azimuth angle 120 • ≤ ϕ ≤ 240 • and elevation angle 10 • ≤ α ≤ 70 • .

•
Step 2: select the reference satellite Using a greedy selection strategy, based on the geometric configuration information provided by the scenario, the PSF of each satellite is estimated by (12), and the reference satellite corresponding to the optimal PSF is selected.

•
Step 3: auxiliary satellite selection The satellite that can be fused with the reference satellite to obtain the best imaging resolution is selected as the auxiliary satellite. Numerical analysis shows that the minimum area of 2-D resolution is to be observed when the two PSF orientations η are nearly orthogonal. Therefore, the angle ψ = η R − η m between the PSF orientation η R of the reference satellite and the PSF orientation η m of the subset of navigation satellites can be calculated, and when η R and η m are orthogonal to each other (ψ = 90 • ), the system achieved the optimal 2-D resolution; when η R and η m parallel (ψ = 0 • or ψ = 180 • ), the 2-D resolution unit of the reference satellite overlaps the 2-D resolution unit of the fusion satellite, the system achieved the worst resolution. Therefore, find the auxiliary satellite in the subset of navigation satellites by an iterative search method, when ψ is closest to 90 • , the corresponding satellite is selected as the auxiliary satellite.
The pseudo-codes of satellite selection strategy are listed as follows:

Algorithm 1: Satellite Selection Strategy
Load the raw data; Calculate receiver position and visible satellites n; for visible satellites n Calculate the azimuth angle ϕ and elevation angle α; Perform threshold comparison for ϕ and α; Obtain a subset of imaging satellites m; end Initiate the processing parameters; for imaging satellites m Calculate the area S m and direction of PSF η m ; Search for minimum area S m ; Search for the best angle ψ; Obtain reference satellite and auxiliary satellite; end

Imaging Algorithm
For better imaging results, we chose the traditional time-domain imaging algorithm, the block diagram of which is shown in Figure 8, Algorithm 2 gives the specific implementation of the pseudo code. The proposed algorithm is an improvement of the (Back projection) BP algorithm, a specific process can be divided into four parts.

1.
Signal Synchronization GNSS-based Bistatic SAR is mainly faced with three synchronization problems: (1) space synchronization; (2) time synchronization; (3) frequency and phase synchronization. The article [35] gives a detailed description of the synchronization problem, and we will not repeat it here.

Range compression and range migration correction
It is well known that the echo channel of the GNSS-based Bistatic SAR system can receive echo signals from multiple satellites. When processing the echo signals, the echo signals from different satellites are superimposed on each other. The pseudo-random code (PRN) in the navigation signal has particularly excellent auto-correlation and cross-correlation characteristics (taking GPS C/A code as an example, the ratio of the auto-correlation peak to the maximum cross-correlation peak is about 24 db). Therefore, the corresponding matched filters can be constructed through the PRN codes of different satellites to achieve range compression without aliasing. First, the carrier phase and code phase information of the satellite's direct signal reconstruct the direct signal, and perform synchronization processing on the reconstruct direct signal and the echo signal, respectively perform range FFT on the reference signal and the echo signal, and convert the signal to time domain and multiplied by the complex conjugate of the reference signal. The result uses Inverse Fast Fourier Transform (IFFT) to complete the range compression of GNSS-bistatic SAR, and the signal is expressed in the following form.
where R x [·] denotes the cross-correlation function between the received and reference signals in the range-time dimension. ω(u − u ta ) is the envelope of the received signal in the azimuth-time dimension, and f c is the carrier frequency of the transmitted signal, and The instantaneous navigation satellite-to-target range is denoted by R T (t), R R (t) is the instantaneous receiver-to-target range, and R B (t) is the instantaneous satellite-to-receiver range. Due to the different equivalent squint angle caused by the linear component of the residual satellite range, the Range Cell Migration (RCM) of the data in the same range cell has a different slope in the azimuth time domain. Although in GNSS-based Bistatic SAR, the transmitted signal is no longer a chirp signal, the signal still has chirp characteristics, so we can use the chirp scaling algorithm to process RCM Correction (RCMC). The filter of RCMC can be expressed as: where f r is the range frequency, transform the time-domain range compression result F(t, u) into the range-frequency domain F( f r , t) through the Fast Fourier Transform (FFT) of the range, multiply it with the range migration correction filter, and then perform Inverse Fast Fourier Transform (IFFT) to complete the range migration correction. The signal is converted back to the time domain, it can be expressed as

Azimuth compression and compensate the Doppler phase
Unlike the RD algorithm, the azimuth compression of the BP algorithm is completed in the time domain, and the specific implementation method is the coherent accumulation of the range compression signal and compensate the Doppler phase in the synthetic aperture time. The azimuth compressed signal can be expressed as where n is the azimuth sampling point.

Fusion imaging results
Using the proposed imaging algorithm, the imaging results of the reference satellite and the auxiliary satellite are obtained respectively. After image registration and amplitude equalization, the final dual-satellite fusion image can be obtained by the noncoherent addition.
The pseudo-codes of fusion imaging algorithm are listed as follows:

Simulation and Discussion
In this section, we demonstrate the validity and feasibility of the proposed method. Potentially, the proposed method can be applied to all constellations of the GNSS system. The 24 h real GPS-L1 orbit data were used for point targets imaging simulation experiment. The precision ephemeris is provided by the International GNSS Service (IGS), which offers three types of precision ephemeris:

IGS Final Precision Ephemeris (IGS Final), the IGS Rapid Precision Ephemeris (IGS Rapid), and IGS
Ultra-Rapid Precision Ephemeris (IGS Ultra-Rapid). In this paper, we chose the IGS Ultra-Rapid to develop the simulation because it has good real-time performance, with an orbital position error of less than 10 cm and a satellite clock error of less than 5 ns. The separation time interval of IGS Ultra-Rapid data is 15 min, it is too long for the sampling interval 1 ms of GNSS-based Bistatic SAR system, Lagrange interpolation was used to complete the data fitting. Table 3 lists the specific parameters, and the distribution of the simulation scene as show in Figure 9. In Figure 10, the corresponding value of the diamond represented the area of the optimal resolution unit obtained by the proposed method, and the corresponding value of the circle represented the area of the global optimal resolution unit obtained by the traversal fusion method. In 24 data sets, the optimal resolutions generated by dual-satellite fusion were significantly better than the optimal resolutions generated by single-satellite. Compared to global optimal dual-satellite fusion resolution, our proposed method generated 15 times the global optimal resolution and 8 times the suboptimal resolution, which showed that the validity and feasibility of the method. In Figure 10, the corresponding value of the diamond represented the area of the optimal resolution unit obtained by the proposed method, and the corresponding value of the circle represented the area of the global optimal resolution unit obtained by the traversal fusion method. In 24 data sets, the optimal resolutions generated by dual-satellite fusion were significantly better than the optimal resolutions generated by single-satellite. Compared to global optimal dual-satellite fusion resolution, our proposed method generated 15 times the global optimal resolution and 8 times the suboptimal resolution, which showed that the validity and feasibility of the method. To show the performance of the proposed method, the area of the optimal resolution units generated by the proposed method were compared with the area of the optimal resolution units of single-satellite. A set of data was selected from 24 data sets for detailed analysis, the data collection time was at 18:00. Using the proposed satellite selection strategy, in Figure 11, there were 4 navigation satellites in the orbit that can be used for imaging, the satellite PRN21 was selected as the reference satellite, and PRN10 was selected as the auxiliary satellite. The single-satellite imaging results of PRN21 and PRN 10 are demonstrated in Figure 12a,b, the dual-satellite fusion imaging result with  n Figure 10, the corresponding value of the diamond represented the area of the op tion unit obtained by the proposed method, and the corresponding value of the c sented the area of the global optimal resolution unit obtained by the traversal fusion metho ta sets, the optimal resolutions generated by dual-satellite fusion were significantly better timal resolutions generated by single-satellite. Compared to global optimal dual-satellite fu tion, our proposed method generated 15 times the global optimal resolution and 8 time timal resolution, which showed that the validity and feasibility of the method. o show the performance of the proposed method, the area of the optimal resolution ated by the proposed method were compared with the area of the optimal resolution un -satellite. A set of data was selected from 24 data sets for detailed analysis, the data colle as at 18:00. Using the proposed satellite selection strategy, in Figure 11, there were 4 navig ites in the orbit that can be used for imaging, the satellite PRN21 was selected as the refer ite, and PRN10 was selected as the auxiliary satellite. The single-satellite imaging resu 1 and PRN 10 are demonstrated in Figure 12a,b, the dual-satellite fusion imaging result To show the performance of the proposed method, the area of the optimal resolution units generated by the proposed method were compared with the area of the optimal resolution units of single-satellite. A set of data was selected from 24 data sets for detailed analysis, the data collection time was at 18:00. Using the proposed satellite selection strategy, in Figure 11, there were 4 navigation satellites in the orbit that can be used for imaging, the satellite PRN21 was selected as the reference satellite, and PRN10 was selected as the auxiliary satellite. The single-satellite imaging results of PRN21 and PRN 10 are demonstrated in Figure 12a,b, the dual-satellite fusion imaging result with our proposed algorithm are shown in Figure 12c, it can be seen that significantly the resolution units of dual-satellite fusion imaging results are smaller than single-satellite.
Sensors 2020, 20, x FOR PEER REVIEW 15 of 21 our proposed algorithm are shown in Figure 12c, it can be seen that significantly the resolution units of dual-satellite fusion imaging results are smaller than single-satellite.  The 2-D resolution was characterized by the area of an approximate ellipse in the 3 dB plane. To better visualize the performance gains resulting from the proposed fusion method, the imaging results have projected onto the 3 dB plane. Figure 13c shows the 3 dB image response obtained the fusion of PRN21 and PRN 10, which area is 853 m 2 , compared with the PRN 21 (1579 m 2 ) in Figure  13a, the area is reduced by 45%. The comparative analysis results fully prove that the proposed method has the excellent resolution improvement performance. our proposed algorithm are shown in Figure 12c, it can be seen that significantly the resolution units of dual-satellite fusion imaging results are smaller than single-satellite.  The 2-D resolution was characterized by the area of an approximate ellipse in the 3 dB plane. To better visualize the performance gains resulting from the proposed fusion method, the imaging results have projected onto the 3 dB plane. Figure 13c shows the 3 dB image response obtained the fusion of PRN21 and PRN 10, which area is 853 m 2 , compared with the PRN 21 (1579 m 2 ) in Figure  13a, the area is reduced by 45%. The comparative analysis results fully prove that the proposed method has the excellent resolution improvement performance.  The 2-D resolution was characterized by the area of an approximate ellipse in the 3 dB plane. To better visualize the performance gains resulting from the proposed fusion method, the imaging results have projected onto the 3 dB plane. Figure 13c shows the 3 dB image response obtained the fusion of PRN21 and PRN 10, which area is 853 m 2 , compared with the PRN 21 (1579 m 2 ) in Figure 13a, the area is reduced by 45%. The comparative analysis results fully prove that the proposed method has the excellent resolution improvement performance.
Sensors 2020, 20, x FOR PEER REVIEW 15 of 21 our proposed algorithm are shown in Figure 12c, it can be seen that significantly the resolution units of dual-satellite fusion imaging results are smaller than single-satellite.  The 2-D resolution was characterized by the area of an approximate ellipse in the 3 dB plane. To better visualize the performance gains resulting from the proposed fusion method, the imaging results have projected onto the 3 dB plane. Figure 13c shows the 3 dB image response obtained the fusion of PRN21 and PRN 10, which area is 853 m 2 , compared with the PRN 21 (1579 m 2 ) in Figure  13a, the area is reduced by 45%. The comparative analysis results fully prove that the proposed method has the excellent resolution improvement performance.

Experimental Verification and Analysis
To demonstrate the validity and feasibility of our proposed fusion imaging formation method, real scene experiments were carried out. To obtain better resolution performance, according to the known accurate ephemeris data, accurate satellite orbit information was obtained, and the following experiment was designed. Figure 14 shows the relevant configuration of the experiment, and Table 4 lists the detailed parameters. The experimental hardware is designed by Beihang University and contains four navigation signal receiving channels. It is an X86 architecture data acquisition and processing platform that supports GPS L1/L5 and BeiDou B1/B3 data acquisition and processing. The right-handed circularly polarized (RHCP) omnidirectional antenna is GPS-702-GG (Novatel, Calgary, AB, CAN) to receive the direct signal of GNSS, which is used to realize system positioning and obtain accurate carrier phase and code phase information to provide reference information for signal synchronization. The echo channel uses a high-gain left-handed circular polarized (LHCP) antenna receives the echo signal from the target area. The LHCP antenna is a 4-array microstrip antenna designed by Beihang University, the specific parameters are shown in Table 4. The data collection location is the roof of Block E of the new main building of Beihang University. The target area is mainly the gymnasium of Beihang University, which is an all-metal roof with a large RCS, that is beneficial for imaging.

Experimental Verification and Analysis
To demonstrate the validity and feasibility of our proposed fusion imaging formation method, real scene experiments were carried out. To obtain better resolution performance, according to the known accurate ephemeris data, accurate satellite orbit information was obtained, and the following experiment was designed. Figure 14 shows the relevant configuration of the experiment, and Table 4 lists the detailed parameters. The experimental hardware is designed by Beihang University and contains four navigation signal receiving channels. It is an X86 architecture data acquisition and processing platform that supports GPS L1/L5 and BeiDou B1/B3 data acquisition and processing. The right-handed circularly polarized (RHCP) omnidirectional antenna is GPS-702-GG (Novatel, Calgary, AB, CAN) to receive the direct signal of GNSS, which is used to realize system positioning and obtain accurate carrier phase and code phase information to provide reference information for signal synchronization. The echo channel uses a high-gain left-handed circular polarized (LHCP) antenna receives the echo signal from the target area. The LHCP antenna is a 4-array microstrip antenna designed by Beihang University, the specific parameters are shown in Table 4.The data collection location is the roof of Block E of the new main building of Beihang University. The target area is mainly the gymnasium of Beihang University, which is an all-metal roof with a large RCS, that is beneficial for imaging.    The experiment was conducted at 20:00 on 20 August 2020. Based on our proposed satellite selection strategy, there are 3 navigation satellites in this area that can be used for imaging (PRN12, PRN15, PRN24), GPS satellite PRN15 was selected as the reference satellite, and PRN24 was selected as the auxiliary satellite. The imaging results of the navigation satellites obtained through the improved BP algorithm are shown in Figure 15a-c. Different images of the gymnasium were obtained under different geometric configurations. The gymnasium is a 100 × 100 m building, subject to the resolution ability of GPS L1 signals and geometric configuration, the obtained image can only reflect that the stadium is a point target in the range direction, and it can be divided into discrete point targets in the azimuth direction. As shown in Figure 15d, the experimental results are consistent with the theoretical calculation. The experiment was conducted at 20:00 on 20 August 2020. Based on our proposed satellite selection strategy, there are 3 navigation satellites in this area that can be used for imaging (PRN12, PRN15, PRN24), GPS satellite PRN15 was selected as the reference satellite, and PRN24 was selected as the auxiliary satellite. The imaging results of the navigation satellites obtained through the improved BP algorithm are shown in Figure 15a-c. Different images of the gymnasium were obtained under different geometric configurations. The gymnasium is a 100 × 100 m building, subject to the resolution ability of GPS L1 signals and geometric configuration, the obtained image can only reflect that the stadium is a point target in the range direction, and it can be divided into discrete point targets in the azimuth direction. As shown in Figure 15d, the experimental results are consistent with the theoretical calculation.   Figure 16 shows the imaging results of dual-satellite fusion, which were generated by the non-coherent addition of the imaging results of the single-satellite. When compared with the area of the minimum resolution unit of Figure 15a, which were reduced about 70.1%, 64.2% and 41.4% in Figure 16a-c, the improvement of the resolution brought dual-satellite fusion was obvious. Compared to the resolution of Figure 16b,c, the Figure 16a was improved most obviously, which was generated by our proposed method. The specific parameters are shown in Table 5, which proves the feasibility and effectiveness of this method. The obtained image is compared with the optical image of the target area. It can be seen from Figure 15d that the stadium can only reflect discrete point targets distributed along with the azimuth direction, due to limited by the range resolution. As shown in Figure 16d, it can be seen that the discrete surface targets revealed in the radar image of the stadium more accurately show the feature information of the stadium, due to the increase of resolution.
can be seen that the discrete surface targets revealed in the radar image of the stadium more accurately show the feature information of the stadium, due to the increase of resolution.  In the experiment, m = 3, the total calculation of the proposed method can be approximated as  ∆S are the experimental measurement that extracts the highlighted partial resolution units from the real image and calculates the area of its 3dB plane approximating the ellipse by Monte Carlo method.
In the experiment, m = 3, the total calculation of the proposed method can be approximated as N a N r log 2 N r + N a N r + N kel N x N y N a * 2, the total computational load for the iterative fusion method can be represented by N a N r log 2 N r + N a N r + N kel N x N y N a * 3, where N a and N r are the number of range and azimuth samples, respectively. N x and N y are the numbers of samples of the imaging scene, and N kel is the length of the interpolation kernel. The experiment is carried out in MATLAB 2018b (MathWorks, Natick, MA, USA) on a computer with AMD R7-3800X 3.9 GHz processors and 64 GB RAM. In the experiment, the time consumed using the proposed algorithm was 1970 s, while that of the iterative fusion method was 2955 s. The experiment demonstrated the feasibility of the proposed method by showing that it achieves good imaging resolution while reducing the computational effort.

Conclusions
In this paper, we focused on the influence of geometric configuration on the fusion resolution of GNSS-based bistatic SAR system. The theoretical analysis and numerical simulation of MPSF are carried out. The influence of geometric configuration and the fusion satellite number on the MPSF are evaluated in combination with point target simulation. It is concluded that dual-satellite fusion is a cost-effective combination. An improved dual-satellite fusion imaging method is proposed. It consists of two parts, a satellite selection strategy and an improved BP algorithm. As demonstrated by simulation and experiment results, the proposed method performed well in terms of both imaging resolution and efficiency. It is worth noting that the method is more effective at long integration times and in airborne mode, where the azimuthal resolution is about 1 m and the area of the fused resolution unit can be reduced to about 1 m 2 . The next step will be a study of airborne mode or multi-system signal fusion with long synthetic aperture time to assess the full potential and challenges of the proposed approach.
In terms of the proposed method itself, we found three application difficulties, which can be improved in the future work.
(1) Multi-satellite fusion has resulted in improved resolution, however, it has introduced a series of side-lope. This will result in a degradation of image quality, however, we plan to use the clean algorithm [29] to solve this problem. (2) The proposed satellite selection strategy is based on the greedy algorithm, and it is easy to fall into local optimum. Therefore, we can optimize the satellite selection strategy to improve the imaging resolution in the follow-up work. (3) The fusion imaging algorithm is realized by an improved BP algorithm, however, it has a large amount of calculation, and we can use an frequency domain imaging method to improve the computational efficiency.

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