Enhanced Passive GNSS-Based Radar Imaging Based on Coherent Integrated Multi-Satellite Signals

Weak reflected signal is one of the main problems in a recent developing remote sensing tool—passive GNSS-based radar (GNSS radar). To address this issue, an enhanced GNSS radar imaging scheme on the basis of coherently integrating multiple satellites is proposed. In the proposed scheme, to avoid direct signal interference at surveillance antenna, the satellites that used as transmission of opportunity are in backscattering geometry model. To coherently accumulate echo signal magnitudes of the scene center in the targeted sensing region illuminated by the selected satellites, after performing the paralleled range compressions, a coordinates alignment operator is performed to the respective range domains, in which, pseudorandom noise (PRN) code phases are aligned. Thereafter, the coordinates aligned range compressed signals of the selected satellites are coherently integrated along azimuth domain so that imaging gain is improved and azimuth processing can be accomplished in only one state operation. The theoretical analysis and field proof-of-concept experimental results indicate that compared to both conventional bistatic imaging scheme and the state-of-the-art multi-image fusion scheme, the proposed scheme can provide a higher imaging gain; compared to the state-of-the-art multi-image fusion scheme, the proposed scheme has a less computational complexity and faster algorithm speed.


Introduction
Global Navigation Satellite System (GNSS) is a satellite system that produces and transmits radio signals for navigation and positioning purposes at global coverage level [1,2]. Up to the year 2020, modernization restructuring of the main GNSS systems, i.e., GPS, GLONASS, Galileo, and Compass (Beidou), will be completed. Throughout existing research works in GNSS area, there already has a wide range of surveys with respect to improving navigation and positioning accuracy on the basis of direct signals. However, in the last decade, the utility of multi-path GNSS signal, known as GNSS reflected signal, gained much attention. Passive GNSS-based radar (GNSS radar) [3,4] is a typical system that uses the reflected GNSS signal as source of opportunity for environmental surveillance. Compared to traditional active radar [5], as there is no need to construct a radar transmission platform, GNSS radar has a lower cost budget and is more flexible for installation under many environmental sensing scenarios. Compared to other passive radar, such as DVT-B based radar [6], because GNSS signals are global coverage and the transmission never failed, GNSS radar can perform all day all weather surveillance. scheme can provide a higher imaging gain and lower computational complexity compared to the state-of-the-art multi-image fusion scheme. This paper is organized as follows. The considered geometry and signal model is given in Section 2. Section 3 analyzed imaging gain and computational complexity of the conventional bistatic imaging scheme and state-of-the-art multi-image fusion scheme, whereas the respective analysis for the proposed imaging scheme is provided in Section 4. Field experimental confirmation of the proposed imaging scheme is shown in Section 5. Section 6 discusses the associated problems in this research and the respective future improvements. Section 7 concludes the whole paper.

The Considered Geometry and Signal Model
In this research, two separated antennas, i.e., direct antenna and surveillance antenna, are used for collecting direct and reflected GNSS signals, respectively. To avoid direct signal interference at surveillance channel for GNSS radar imaging, similar to [29], the satellites in the geometric position of backscattering are chosen as sources of opportunity. The geometry model is shown in Figure 1. In Figure 1, for the ease of calibrating range migration when receiver is moving within certain trace for forming SAR image, both direct and surveillance antennas are mounted on the same platform. The signals received at direct antenna are used for synchronization, whereas the signals collected by surveillance antenna are used for radar imaging. At GNSS receiver, the received signals are digitized, downconverted to base-band and formed into range and azimuth domains. The respective mathematical model with respect to received direct signal for each satellite is given as where i represents the index for each satellite; A i d represents the signal magnitude; C presents PRN code; D represents navigation message; t represents range domain, which is upper bounded by the length of PRN code; u represents azimuth domain, which is upper bounded by data collection duration or receiver moving duration; τ i presents the transmission delay between receiver and each satellite; ω i d represents Doppler frequency of each satellite; φ i d represents carrier phase of each satellite; and n d represents noise at direct antenna. The parameters ω i d and φ i d can be considered as constants in the same range domain.
Each reflected signal can be considered as delayed version of the respective direct signal, which can be expressed as Equation (2) where k represents the index of each reflected signal at range domain, τ k r i represents the reflected signal delay compared to the respective direct signal each satellite, ω i r represents Doppler frequency of reflected signal, φ i r represents carrier phase of reflected signal, and n r represents background noise at surveillance antenna. For stationary target, ω i r = ω i d , whereas for moving target, ω i r − ω i d represents Doppler frequency caused by the object velocity. The parameters ω i r and φ i r can be considered as constants with the same range domain as well.

Analysis of Conventional Bi-Static Imaging Scheme and State of Art Multi-Images Fusion Scheme
GNSS radar is a passive radar, the transmitter and receiver should be located on separated platforms, thus only bistatic imaging scheme and multi-static imaging scheme are appropriate for such kind of system. In the GNSS radar-related research, bistatic imaging scheme [8] is regarded as a conventional scheme, which consists of the stages signal synchronization and radar imaging for individual satellite. Multi-static imaging scheme [9,13,16,18,19], known as the multi-image fusion scheme, is the state-of-the-art scheme, which is functions primarily on the basis of fusing multiple bistatic GNSS radar images. The detailed analysis of these two schemes are given as follows.
First, direct signal synchronization is carried out by tracking the received direct signal (1) from all the visible satellites. Using decoded navigation message, the satellite that satisfy the backscattering geometric model as Figure 1 is selected. In bistatic imaging, only one satellite that satisfied the geometric model in Figure 1 is selected as transmission of opportunity. Based on tracked code delay, carrier phase, and navigation bits for the selected satellite, the local replica is generated as Thereafter, for imaging stage, range compression is performed by correlating s i r with s i m per range domain along azimuth. The range compressed signal can be expressed as where l r denotes the range samples used for compression, N s denotes the respective samples quantity, * represents conjugate operator, represents the bistatic range compressed carrier phase of the i-th satellite. The value 2N s in R i c is due to the fact that the quantity of samples is doubled after performing correlation operator.
For stationary object imaging, generally SAR image is formed. For GNSS-SAR image formation, azimuth matched filter is obtained on the basis of the result in Equation (4) along azimuth domain u. Azimuth compression is identically azimuth matched filtering, which can be expressed as where l u represents azimuth samples for compression and M u denotes respective azimuth samples quantity. The value 2M u in T i S occurs because of the impact of azimuth correlator as well. For moving object, normally RD map is generated. The respective azimuth processing is carried out by performing FT of Equation (4) along azimuth domain, which can be expressed as where ω represents azimuth frequency caused by the movement of object. To obtain GNSS-SAR image or GNSS radar RD map, coordinates in both range and azimuth domains are transformed into distance domains, and absolute operator |·| is applied on Equation (5) or Equation (6).
As for the state-of-the-art multi-image fusion scheme, more than one satellite are used for imaging. First, multiple bistatic images on the basis of each satellite are generated. Then, to transform the coordinates into distance domains, pseudo-range and elevating angle of each satellites are obtained on the basis of decoded navigation message. Thereafter, the scaling factor for image alignment is generated on the basis of the coordinates transformed results, which can be seen in detail in [19]. Assuming the scaling factor for performing alignment for each satellite is ζ i , the generation of GNSS-SAR image and GNSS radar RD map on the basis of multi-image fusion scheme can be expressed as Equations (7) and (8), respectively, where m represents the quantity of satellites used for multi-image fusion. We investigate imaging gain for conventional bistatic imaging and multi-image fusion schemes. In bistatic GNSS-SAR imaging scheme, after performing Equation (4), the range compression gain for individual satellite can be derived as G rc = N s ; For azimuth compression Equation (5), the gain in the total dwell time can be derived as G total = (N s ) 2 · M u ; thus, the gain for azimuth compression can be expressed as G ac = (N s ) 2 ·M u N s = N s · M u . Therefore, the imaging gain for bistatic GNSS-SAR imaging scheme can be expressed as which is the same as G total . For bistatic GNSS radar RD map formation, the gain in range compression is the same as bistatic GNSS-SAR imaging scheme. In terms of azimuth FT, the respective gain can be derived as G a f = M u . Thus, the imaging gain can be expressed as In terms of the multi-image fusion scheme, it can be regarded as the non-coherent integration of Equation (5) or (6) with coordinates alignment. Therefore, it can be easily derived that for multi GNSS-SAR images fusion, the gain is For multi GNSS radar RD maps fusion, the gain is From Equation (9) to (12), it can be seen that compared to bistatic imaging scheme, the gained strength for multi-image fusion is √ m larger. Computational complexity in this paper is studied on the basis of number of operations. In signal synchronization stage of both bistatic imaging scheme and multi-image fusion scheme, as all the satellites are used for processing, the complexities are the same. Thus, only the complexity during imaging stage is considered. In terms of bistatic imaging, for local replica generation, the number of operations in code modulation and carrier modulation is the same as N s × M u , thus the respective complexity is derived as O (2N s × M u ). For range compression state, there exists number of multiplications N 2 s × M u and number of additions N s × M u , thus the respective complexity is O (N s (N s + 1) × M u ). For azimuth compression for stationary object imaging, as the quantity of samples at range domain is doubled after performing compression, the number of multiplications is derived as M 2 u × 2N s and the number of additions is derived as M u × 2N s ; thus, the respective complexity is derived as O (M u (M u + 1) × 2N s ). Thereafter, the accumulated complexity for bistatic GNSS-SAR imaging stage can be derived as As for bistatic RD map generation for moving object detection, the complexity for local replica generation and range compression is the same as bistatic GNSS-SAR imaging. However, in the azimuth FT state, there exists number of multiplications 2N s × M u 2 log 2 M u and number of additions 2N s × M u log 2 M u . Thus, the respective complexity is derived as O (3N s × M u log 2 M u ). Therefore, the accumulated complexity for RD map generation stage can be derived as For multi-image fusion scheme, through a similar theoretical analysis, the number of operations for local replica generation, range compression, azimuth compression for stationary object detection, and azimuth FT for moving object indication can be derived as respectively. The number of operators for both coordinates alignment and images combination with respect to GNSS-SAR imaging can be derived as 2N s × 2M u × m, whereas for RD map forming, the respective number of operations can be derived as 2N s × M u × m. Therefore, the accumulated complexity for multi GNSS-SAR images fusion can be derived as whereas for multi GNSS radar RD maps fusion, the complexity can be derived as

The Proposed Imaging Scheme Using Coherent Integrated GNSS Signals
To improve imaging gain and reduce computational burden than the state-of-the-art multi-image fusion scheme, a new imaging scheme on the basis of coherently integrating multiple GNSS satellites signals is proposed. The main principle of the proposed scheme is to coherently integrate coordinates aligned multi-satellites range compressed signals based on the compressed carrier phase difference. Thereafter, azimuth compression or azimuth FT can be performed by only once-through operation. The detailed analysis can be seen as follows.
The signal synchronization stage and range compression stage for the proposed imaging scheme is the same as multi-image fusion scheme. The criterion for selecting the multiple satellites as sources of opportunity should satisfy the model in Figure 1 as well. The initial range compressed signal without coordinates alignment is the same as Equation (4).
Thereafter, range coordinates alignment is performed. As it is hard to obtain the spatial information of the glisten region on the passive radar platform without generating full GNSS-SAR image or full GNSS radar RD map, unlike [19], the range coordinates alignment in the proposed scheme is carried out based on the synchronized carrier phase difference of direct signals among the selected satellites. The detailed derivation can be seen as follows. • At first, bistatic range distance difference of the selected satellites for imaging is considered. We use one of the selected satellites as a benchmark and mark it as 0-th satellite. Assume, in the 0-th satellite, the distance between satellite and object is R 0 t and between object and receiver is R 0 r ; then, the respective bistatic range can be calculated as d 0 = R 0 r + R 0 t . Assume for the satellite i, the distance between satellite and object is R i t , and between object and receiver is R i r ; similarly, the respective bistatic range is calculated as Because the distance between object and receiver is not related to satellite position, we can have that R i r = R 0 r . Thus, the bistatic range difference is derived as We investigate the impact in the changes of range distance between object and receiver with respect to GNSS satellite elevation angle. For the ease of analysis, we plot a respective schematic diagram in Figure 2. In Figure 2, R b represents the distance between satellite and receiver; h represents vertical distance between satellite and earth surface; θ d and θ r represent elevating angles at receiver and at object, respectively; and ∆θ represents the difference between θ d and θ r . Based on Figure 2, the changes of range distance between object and receiver R r with respect to θ r and ∆θ is calculated as We study the changes of R r when ∆θ reaches its minimum value 1 • . Under the circumstance, GPS satellite is used as an example, in which the average vertical distance h = 22,200 km, the relationship between R r and θ r is simulated in Figure 3, where the general interval of elevating angle for GNSS satellite is between 10 • and 80 • [31]. From Figure 3, note that to reach the level of elevation angle change by only 1 • , the change of R r should be at the level more than 10 7 meter. As GNSS radar is a passive radar, it is unlikely that the sensing range can reach that level. Therefore, under majority circumstance, it can be regarded that θ d = θ r . Then, the expression Equation (17) can be transformed into Equation (19) as follows.
where R i b and R 0 b represent the base-line distance between receiver and the 0-th and the i-th satellite, respectively; c represents signal transmission speed; and f represents signal carrier frequency. The carrier phase values φ i d and φ 0 d of the i-th and 0-th satellites can be obtained from tracked results of Equation (1).
• Second, the necessity for performing coordinates alignment is evaluated. If ∆d i < c N s × 10 −3 , where c N s × 10 −3 represents range distance between two sampling point within the period 1 ms of the PRN code for civil purpose, there is no need to perform coordinates transform; otherwise, perform coordinates alignment.
• Third, the aligned coordinates using the 0-th satellite as a benchmark can be derived as Equation (20) After coordinates of the satellites are well aligned with the 0-th satellite, range compressed signals with respect to the satellites as illuminator of opportunity are accumulated coherently along azimuth, which can be expressed as After performing Equation (21), the scene center of the point spread function (PSF) with respect to illuminated region can be coherently accumulated. However, the illuminated ambiguity region of PSF with respect different satellites will be different. This will negatively impact on range resolution after performing Equation (21). As the main aim in this paper is the preliminary feasibility investigation with respect to imaging on the basis of coherently integrated multiple GNSS satellites, the range resolution problem is not specifically concentrated. For stationary object imaging, azimuth compression is carried out based on the result Equation (21) along azimuth domain, which can be expressed as For moving target detection, on the basis of Equation (21), azimuth FT is performed as follows.
Applying absolute operator on Equation (22) or (23), final GNSS-SAR image or GNSS radar RD map with respect to coherently integrated satellites can be obtained.
In summary, the work-flow of the proposed imaging scheme is given as Algorithm 1.

Algorithm 1
The work-flow of the proposed imaging scheme 1.
Performing signal synchronization, selecting the GNSS satellites that satisfy the model as Figure 1.

2.
On the basis of the selected satellites, generating local replica for each satellite as Equation (3).

3.
Performing range compression for each selected satellite independently as Equation (4).

4.
Using one of the selected satellite as a benchmark, mark it as 0-th satellite, extracting carrier phases of direct signals on the basis of the selected satellites.
, performing coordinates alignment, otherwise directly jump to step 6.

5.
Obtaining the coordinates alignment factor as Equation (19), and performing range coordinates alignment along azimuth as Equation (20).

6.
Coherently combining the range compressed signals along azimuth direction of the selected satellites as Equation (21).
Imaging gain and computational complexity of the proposed imaging scheme are analyzed. The gained strength for each satellite in range correlation operation for the compression is the same as the bistatic imaging scheme. However, after performing Equation (21), the gain can be derived as G r_proposed = N s · m, after performing Equation (22), and the azimuth gain can be derived as . Thus, the gain for GNSS-SAR imaging under the proposed scheme can be derived as After performing Equation (23), the gain is derived as G a f _proposed = M u · m. Thus, the imaging gain for RD map generation under the proposed scheme can be derived as Comparing Equation (24) with Equation (11), it can be seen that the gained strength in GNSS-SAR imaging under the proposed scheme is m 3 2 higher than multi-image fusion scheme; comparing Equation (25) with Equation (12), it can be seen that the gained strength for GNSS radar RD map generation under the proposed scheme is m whereas for GNSS radar RD map generation, it is Comparing Equation (26) with Equation (15), it can be seen that the proposed scheme has less complexity by 2N s × M u × (2m + (M u (m − 1) − 2)) number of operations than the multiple GNSS-SAR images fusion scheme; comparing Equation (27) with Equation (16), it can be seen that the proposed scheme has a less complexity by 3N s × M u × (m − 1) log 2 M u than multiple GNSS RD maps fusion scheme.

Proof-of-Concept Field Experimental Confirmation
To validate the proposed scheme for enhancing GNSS radar imaging, field proof-of-concept experiments were carried out on the basis of software defined radio (SDR) GPS C/A code receiver. The same experimental equipment as first author's previous work [32] is employed in this research. The equipments at receiver end are illustrated in Figure 4.   In Figure 4a, the direct antenna is right hand circular polarization (RHCP), which faces the sky to collect direct GPS signal for performing signal synchronization. The surveillance antenna is left hand circular polarization (LHCP), which used for collecting reflected GPS signals for the targeted surveillance region. For maintaining time synchronization, both direct and reflected signals are saved by the same radio frequency (RF) front end with two separated channels. The model of the RF front end is ET09/C, which produced by ip-solution company (http://www.ip-solutions.jp) . The interface of the software for signal data collection is given in Figure 4c. The collected raw signal data were processed on commercial computer platform. The parameter values for the experiments is given in Table 1. where T denotes PRN code period. Based on the experimental set-up and parameters, two representative scenarios, i.e., land imaging scenario with two strong reflectors and ocean target detection scenario, are considered. The detectability of the proposed scheme is examined based on imaging gains. The respective analyses can be seen as follows.

Land Imaging Scenario With Two Strong Reflectors
Under this scenario, GNSS radar worked in SAR imaging mode. The satellites were considered to be the stationary transmitters. The counterclockwise movement of the rotator in Figure 4a is employed for performing azimuth angular trace of surveillance antenna for aperture synthetic. The length of azimuth trace is 60 • , while the moving duration is 1 min. Within the trace, optical image of the surveillance region is shown in Figure 5.
In Figure 5, range distance between GNSS radar receiver and the slope is approximately 30-40 m, whereas the azimuth distance between the two reflector is~30 • . The sizes of the two reflectors are 0.2 m × 0.2 m. According to the decoded navigation message from direct signal synchronization, the satellites GPS PRN 15 and PRN 29 were used as source of opportunity, as they were only the satellites in the geometric position of backscattering under this scenario. First, the authors estimated bistatic range distance difference of the two satellites based on the carrier phase values from synchronization. The carrier phase difference can be seen in Figure 6a, whereas the coordinates transform results in distance domain can be seen in Figure 6b. From Figure 6, it can be seen that the error caused by the position of two satellites is larger than the range distance represented by the distance between two sampling points 18 m. Thus, it is required to perform range domain alignment. As PRN 29 is further than PRN 15, on the GNSS radar platform in this experiment, after performing the separated range compression stages, the authors used the bistatic range distance of PRN 15 as a benchmark, whereas the bistatic range distance of PRN 29 minus the values shown in Figure 6b per range domain along azimuth time for the calibration. Thereafter the range compressed signals of the two satellites are coherently integrated as Equation (21) for performing azimuth compression. The obtained GNSS radar image with respect to the proposed scheme can be seen in Figure 7. For comparison, the images obtained by multi-image fusion scheme and bistatic images of the two satellites are given in Figures 8 and 9, respectively. We investigated imaging gain by examining the highest pixel intensity. From Figures 7-9, it can be seen that as the integration of bistatic results were employed, the imaging gain in both Figures 7 and 8 are higher than Figure 9. Comparing Figure 7 to Figure 8, because the proposed imaging scheme works in the mode that coherently integrating the coordinates aligned range compressed signal, the imaging gain is approximately 2.83 higher than the state-of-the-art scheme multi-image fusion, where the value m in this scenario is 2.

Azimuth domain (mini-seconds)
The computational efficiency for this experiment is studied on the basis of numbers of operations during imaging and the respective algorithm speeds. From Table 1, we have that the number of samples in each range domain is 16368. As two satellites were used as signal of opportunity, the PRN code period is 1 ms and the moving duration for signal collection is 1 min, we can have that the number of samples in azimuth domain is 60 × 10 3 . Therefore, under the proposed imaging scheme, there exists 2 × 16368 × 60 × 10 3 × 2 = 3.93 × 10 9 numbers of operations during local replica generation. For imaging range compression stage, the number of operations is 16368 × 16369 × 60 × 10 3 × 2 = 3.22 × 10 13 . The number of operations during coordinates alignment and range signal coherently accumulation are 3.93 × 10 9 as well. As coordinate aligned range compressed signals are coherently integrated along azimuth, only one azimuth compression stage is needed. From the derivation with the experimental parameter values, the numbers of operations during one azimuth compression stage is 60 × 10 3 × 60 × 10 3 + 1 × 2 × 16368 = 1.18 × 10 14 . Then, we can have that the accumulated number of operations during the procedure is approximated 1.20 × 10 14 . For multi-image fusion scheme, the difference under this scenario is that there exists two azimuth compression stages for the satellite GPS PRN 15 and PRN 29, respectively. Through a similar analysis, the number of operations during the imaging procedure can be derived as 2.02 × 10 14 . For bistatic imaging scheme with either GPS PRN 15 or PRN 29, as only one satellite is employed as transmitter of opportunity, the number of operations will be deducted by half during azimuth compression compared to multi-image fusion scheme, and deducted by half during local replica generation state and range compression state, compared to both the proposed imaging scheme and multi-image fusion scheme as well. Meanwhile, there is no coordinate alignment step. The number of operations with respect to the three schemes during imaging stage are summarized in Table 2.  The algorithm speeds are studied on the basis of machine running time under the same computer environment, which are concluded in Table 3.  Table 2 to Table 3, it can be seen that the proposed scheme has less computational complexity than multi-image fusion scheme. Meanwhile, the proposed scheme is faster than multi-image fusion scheme. However, as both the proposed scheme and multi-images scheme need to generate local replica and perform range compression of the signals from the two satellites independently, their computational complexity and the time spends on the procedure for imaging are higher than bistatic imaging scheme.

Ocean Moving Object Detection Scenario
This subsection considers the ocean moving object surveillance scenario, in which RD maps are forming. The experiment was carried out at Hong Kong Cyberport. The optical image for this scenario is given in Figure 10.
In this scenario, we aim to form RD maps of the ocean ferry. To form backscattering geometric model for imaging, on the basis of decoded navigation message, the satellites GPS PRN 15 and GPS PRN 24 are used as source of opportunity. From the measurement by laser range finder, the distance between the ferry and surveillance antenna is approximately 910-930 m. Both direct and surveillance antennas are fixed by the rotator (see Figure 4) for raw signal data collection. The duration for signal collection is 1 min. Therefore, the signals samples along azimuth domain is 60 × 10 3 as well. Based on the respective signal collection set-up, first, the bistatic range distance difference along azimuth domain between the two satellites are extracted based on synchronized carrier phase value, which can be seen in Figure 11. From Figure 11, it can be seen that the bistatic range distance difference between the two satellites is less than the distance between two sampling points 18 m. Thus, the respective range compressed signals can be directly integrated without performing coordinate alignment. After performing the coherent accumulation of the range compressed signals from the respective satellites, azimuth FT was carried out for obtaining RD maps. The generated RD maps based on the proposed imaging scheme, multi-image fusion scheme and bistatic imaging scheme are given in Figure 12. For the ease of comparison, the results are illustrated in three-dimensional mode. Comparing the highest pixel intensity among each image in Figure 12, we can see that in the proposed imaging scheme, as the range signals are coherently integrated, indeed Figure 12a has the highest imaging gain. Because multi-image fusion scheme is a non-coherent integrated scheme, the imaging gain in Figure 12b is not as high as Figure 12a.
Similar as Section 5.1, the efficiency of this experiment is investigated based on numbers of operations and algorithm speeds extracted from machine running time. Based on the parameter values in Table 1 and the data collection duration, for the proposed imaging scheme, the number of operations with respect to local replica generation is the same as the respective results in Section 5.1. For the ease of comparison, range compression was carried out at frequency domain in this experiment. Through analysis, it has the number of operations 2 × 3 2 × 16,368 ×log 2 (16,368) = 687,380 for range FT and range inverse FT, respectively. The number of operations for range multiplication for the compression is 2× 16,368 = 32,736. There is no operations for coordinates transform, whereas the number of operations for coherently accumulating range compressed of the selected satellites is the same as the respective result in Section 5.1 as well. As range compressed signals of GPS PRN 15 and PRN 24 are accumulated, only one stage azimuth FT is necessary. The number of operations with respect to complex multiplication and addition during performing azimuth FT in each range position are 60×10 3 2 log 2 60 × 10 3 ≈ 4.76 × 10 5 and 60 × 10 3 · log 2 60 × 10 3 ≈ 9.52 × 10 5 , respectively. Therefore the accumulated number of operations during imaging processing for RD map forming is derived as 2559500. For multi-image fusion scheme, there exists two azimuth FT processing for GPS PRN 15 and PRN 24, respectively. Thus, there will exist additional number of operations 1428000. The number of operations with respect to the proposed imaging scheme, multi-image fusion scheme and bistatic RD mapping scheme during the imaging stage are summarized in Table 4. The algorithm speeds based on machine running time are illustrated in Table 5.  Tables 4 and 5, the same verdict as Section 5.1 can be concluded: the proposed imaging scheme has less computational burden and faster than multi-image fusion scheme.

Discussion
In this paper, due to the fact that it mainly focuses on the feasibility testing with respect to the proposed imaging scheme based on integrating satellites coherently, only bistatic range is considered. Thus, only the information carrier phase difference is used for range coordinates alignment. However, to estimate the objects location more precisely, the parameter values that elevating angle and pseudo-range between each satellite and receiver are still required. Therefore, in future, the authors aim to improve range coordinates alignment stage in the proposed imaging scheme for obtaining the objects range position on GNSS radar image more precisely and with a lower computational complexity than the respective state in multi-image fusion scheme.
Meanwhile, the field experiments are carried out on the ground-based GNSS radar. Under many circumstances, only two satellites are satisfied backscattering geometric position as Figure 1. Thus, to maximumly avoid direct signal interference at surveillance channel, only two satellites are used for the coherent accumulation. In the future, the authors aim to perform experiments on airborne GNSS radar platform. As all the satellites are satisfied under the geometric model, as in Figure 1, all of them can be employed as transmitter of opportunity, in which, the advantages of the proposed imaging scheme will be more significant in field implementations. At the same time, as the PSFs around the scene center of the same target illuminated by different satellites are not the same, although the magnitude of scene center can be improved by coherently accumulating range compressed signals, the range resolution will be negatively impacted. Therefore, in future, the authors will further improve the proposed scheme by combining the range resolution enhanced mechanism in the first author's previous work [25]. In this remit, a range resolution enhanced GNSS-SAR or GNSS radar RD map can be obtained on the basis of coherently integrated multi-satellites.
In addition, for this research, the authors only have GPS C/A code SDR receiver as shown in Figure 4 for field experimental testing. In the future, the authors will apply the proposed imaging scheme to GLONASS, Galileo, and Beidou (Compass) signal receivers, in which, the adaptability of the proposed scheme can be enhanced.

Conclusions
A new imaging scheme on the basis of coherently integrating multiple GNSS satellites is proposed in this paper. In the proposed scheme, the satellites that satisfy backscattering model is selected as sources of opportunity. Based on the synchronized carrier phases of the selected satellites, range coordinate alignments among satellites are performed after performing range compressions independently. Thereafter, the coordinate aligned range compressed signals are coherently accumulated along azimuth domain, in which azimuth compression can be completed in only once-through operation. Both theoretical analysis and field proof of concept experiments show that compared to both conventional bistatic imaging scheme and state-of-the-art multi-image fusion scheme, the proposed imaging scheme can have a higher imaging gain; compared to multi-image fusion scheme, the proposed imaging scheme has a lower computational complexity. In conclusion, the proposed scheme will be more suitable for GNSS radar image formation under weak reflected signal circumstance.