Moving Target Detection and Parameter Estimation Using BeiDou GEO Satellites-Based Passive Radar With Short-Time Integration

Owing to the global coverage and high visibility of global navigation satellite systems (GNSS), passive radar using GNSS satellites as transmitters of opportunity has shown great potential for moving target detection, and its possibility to detect short-range maritime targets has been verified. Nevertheless, the existing methods always require tens of seconds for echo integration, which is not ideal for the detection of targets with high maneuverability. To solve this problem, a bistatic configuration using BeiDou geostationary earth orbit (GEO) satellites as illuminators of opportunity, which is suited for short-time coherent integration of echo signals, is presented. In addition, the corresponding signal processing method is introduced. To achieve a more continuous trajectory of the target, a method used to fuse range-Doppler maps obtained individually from multiple BeiDou GEO satellites is also described. The effectiveness and functionality of our proposed algorithms are demonstrated with an experiment in which a ferry was used as a moving target and its real-time location data were adopted as a reference.


I. INTRODUCTION
S INCE global positioning system (GPS) became fully operational in 1993, the idea of using the reflected signals of navigation satellites for remote sensing applications has been a hotspot research topic. Generally speaking, there are two main branches of global navigation satellite system (GNSS) remote sensing used for earth observation, i.e., GNSS reflectometry and GNSS radio occultation. Meanwhile, from the perspective of radar, employing GNSS signals as illuminators of space-surface synthetic aperture radar (SAR) and multistatic passive radar systems is also of great potential [1].
Back in 1988, Hall and Cordey [2] proposed the concept of multistatic scatterometry using GPS signals. In 1996, Komjathy et al. [3] from National Aeronautics and Space Administration Langley Research Center developed the idea of using GNSSreflected signals for earth observation applications.
Since then, numerous works have been performed to implement this technique. From oceanography [4]  monitoring [5], and even for soil moisture monitoring [6], GNSS-reflected signals are widely used. At the same time, the concept of a passive radar system using the GPS signal as the illuminator was proposed in 1995 [7]. The importance of passive radar systems using GNSS signals as illuminators of opportunity is highlighted due to the high coverage and availability of GNSS satellites. Nevertheless, there are still many difficulties to overcome for target detection with GNSS-based passive radar systems. Since GNSS signals are not designed specifically for radar applications, the echo power is always too weak to detect directly. Even worse, when the target to be detected is in motion, the received signals will be more difficult to integrate due to the unpredictable migration of range and Doppler shift. Moreover, the bistatic reflection, which keeps changing during the observation period due to the continuous motion of the satellite and target, can also influence the target's radar cross section along with the power of the echo signal, and further affect the obtained signal-to-noise ratio (SNR) significantly [8], [9]. All the aforementioned drawbacks are prevalent when the GNSS-based passive radar is employed for target detection.
To solve these problems, numerous solutions for moving target detection (MTD) have been proposed. A commonly used method is to integrate the echo signals over the dwell time, and it originates from the MTD technique adopted by conventional radar [10], [11], [12]. Traditional MTD works well when the target moves with a constant radial velocity and does not travel across multiple range cells. The fundamental principle of MTD is to integrate all the received energy from the target into specific range and Doppler cells in the range-Doppler (RD) domain. The index of the focused range cell stands for the distance between the receiver and the target, and the Doppler shift corresponds to the radial velocity of the target. During one coherent integration interval (CPI), if the target moves across multiple range gates, whose size depends on the sampling rate of the intermediate frequency (IF) signal and the signal bandwidth, the received energy from the target would smear to those cells [11]. The above phenomenon is known as range cell migration (RCM). Similarly, the Doppler cell migration (DCM) can be found when the echo energy spread over different frequency cells in the azimuth-frequency domain. In traditional MTD, both RCM and DCM are not compensated, and therefore, if the noncooperative target moves across multiple range or Doppler cells in one CPI, the performance of MTD will deteriorate. Hence, these techniques adopted to eliminate both RCM and DCM should be considered here [13]. So far, various novel algorithms used to correct RCM and DCM have been proposed, and moving targets have been detected successfully at short ranges with a GNSS-based passive radar system using these algorithms. Further applications such as target localization [14], [15], [16], [17], [18] and space-surface bistatic SAR imaging [19] can also be accomplished.
For the compensation of RCM, the most straightforward method is simply removing the phase term related to the radial velocity of the target from the original complex signal. In most cases, however, the velocity is unknown, and the sounding system cannot even determine whether the target exists or not. Hence, it is necessary to compensate for the RCM without prior velocity information. An algorithm used in such a case is the radon-Fourier transform (RFT), which is capable of detecting the target and estimating its motion parameters in the case of translational motion. The possibility of using this method to detect maneuvering targets was examined with simulation results in [20]. A more universal version, referred to as generalized radon-Fourier transform (GRFT), was also developed to be used in these cases, in which the target has higher order motion parameters than the case of simple translational motion. The performance of GRFT has been verified by both simulated and experimental results [21], [22], [23]. However, both RFT and GRFT are essentially a searching strategy to find out the correct motion parameters of the target. Therefore, they are computationally expensive, especially for the target motion expressed by higher-order polynomials. In addition, Keystone transform (KT) is another effective method to eliminate RCM. It can compensate RCM without prior velocity information and reduce the algorithmic complexity of RCM correction, which is its key merit compared with RFT and GRFT. Note that, every possible radial velocity of the target should be verified within a given searching range in RFT and GRFT. KT has been widely used in MTD with conventional pulsed radar but has rarely been applied in GNSS-based passive radar systems [24], [25], [26], [27].
For DCM correction, the problem is more complicated because echo signals need to be transformed into the frequency domain. The DCM is produced when the target moves with a nonzero radial acceleration. To compensate for the DCM, the traditional discrete Fourier transform (DFT) cannot determine the Doppler shift accurately because the signal frequency keeps changing. Numerous methods have been proposed to solve this problem, one of which is known as joint time-frequency analysis. Among these joint time-frequency methods, short-time Fourier transform (STFT) has been widely employed, but its performance is highly dependent on the SNR of the raw echo data and cannot be used in low-SNR scenarios [26], [28], [29], [30]. Wavelet transform (WT) can be used to project the echo signal into a set of basis functions named wavelets, giving out the frequency components and corresponding occurrence time [31], [32]. However, the performance of WT relies on the SNR of each frequency component, so it has similar limitations as STFT. Another method called fractional Fourier transform (FrFT) is more suitable for solving this problem. If a target's motion can be represented by a second-order polynomial, the echo signal in the azimuth-time domain will be approximated by the linear frequency modulation (LFM) signal, and FrFT is designed to be used in this case [33], [34], [35]. For example, FrFT was employed to compensate for DCM in [34], and the simulation and experimental results showed that the ship targets can be detected successfully. However, a CPI of 2.5 s was adopted, which is not short enough for the detection of targets with high maneuverability. In addition, a Dechirping algorithm (DA) also can be applied when the echo signal in the azimuth-time domain can be approximated by the LFM signal. Its principle is similar to FrFT, however, has lower computational complexity due to the more concise and efficient searching strategy for the rate of frequency modulation [36], [19]. In [19], the DA was used to compensate for the in-frame Doppler migration, and then the physical target length was estimated after implementing the inverse synthetic aperture radar technique. Note that, a CPI of 3 s was adopted during their experiments. Therefore, it can be found that a long CPI was usually employed in these previously reported experiments [37], [38], [39], [40], [41]. Long CPI contributes to high SNR naturally but only works when DCM and RCM are corrected properly. However, it is difficult to achieve this task when the target has complicated motion parameters, e.g., the target with high maneuverability. To handle the problem of target motion in this case, a refined segmentation strategy on the recorded echo data should be applied. After this operation, the target motion during the period of each data segment can be simplified. In addition, almost all the experimental results were verified with automatic identification system (AIS) data of the maritime target. The AIS data, however, have an unstable and low update rate, making it can hardly be regarded as a reliable reference.
In this article, a novel method is proposed to integrate the echo signal with a short CPI, and the true movement trajectory of the detected target is used for comparison and verification. In this method, the KT is first implemented to resample data in the range-frequency domain and correct the linear RCM. Then, DA is used to find the chirp rate along the azimuth time domain to compensate for the DCM. After RCM and DCM corrections, the RD map is obtained to show the echo intensity and the radial velocity at the same time. Besides the processing of energy focusing, the radon transform algorithm (RTA) is also adopted to estimate the radial acceleration of the detected target using the bistatic range extracted from the RD map. To verify that our algorithms can effectively detect the moving target and estimate its parameters with short-time integration, the experimental results are compared with the real-time location data, which were calculated from the actual movement trajectory of the target. Finally, an algorithm used to fuse RD maps acquired separately from different BeiDou geostationary earth orbit (GEO) satellites to better present the movement trajectory of the target is also introduced. The remainder of this article is organized as follows. In Section II, the signal model is presented. In Section III, the signal processing method is described in detail. Here, the correction of RCM and DCM is introduced, and methods used to estimate the target's motion parameters are proposed. In Section IV, the experimental results are described. Section V focuses on multisatellite data fusion and Section VI concludes the article.

II. SIGNAL MODEL
The transmitters of opportunity adopted in our study are BeiDou GEO satellites, and the navigation signal broadcast in the B3 in-phase (I) band is used as the illuminating signal. According to the official interface control document (ICD) of the BeiDou navigation satellite system (BDS) [42], the BeiDou B3I components can be defined as follows: (1) where f 0 is the carrier frequency, ϕ B3I is the initial phase, A B3I represents the amplitude, D B3I is the data code, and C B3I is the ranging code.
It is well known that the polarization of BDS transmitted signal is right-handed circular polarization (RHCP). If the incident angle of the BDS signal is larger than Brewster's angle of the scattering medium, then the polarization of the BDS signal will be changed to left-handed circular polarization (LHCP). Since the observed target has metal surfaces, the LHCP antenna can receive stronger reflected signals from the target. Meanwhile, the LCHP antenna can also suppress the RCHP direct signal. Therefore, in our bistatic radar system, two antennas were employed. One is an RHCP antenna connecting to a heterodyne channel (HC) to receive the direct signal, and another is an LHCP antenna connecting to a radar channel (RC) to receive the echo signal. In the HC, the received signal can be expressed as where τ d denotes the delay of the direct signal. The echo signal scattered by a point target can be written as where τ r is the delay of the reflected signal. For simplicity, the envelope of the S B3I_d (t) or S B3I_r (t) is denoted as S(t), i.e., S (t) = A B3I D B3I C B3I (t). In the multi-channel receiver, the received signal is downconverted to the IF and demodulated to a complex signal. For direct signal S B3I_d (t), it is separately mixed by cos[2πf m (t)] in the I channel and by sin[2πf After downconversion and demodulation, the received signal in HC can be described as Similarly, the received signal in RC after down-conversion and demodulation can be described as Here, a new parameter τ delta is introduced to represent the delay difference between the reflected and direct signal and can be written as where R st is the range between the satellite and target, R tr represents the range between the target and receiver, R sr is the range between the satellite and receiver, and R delta represents the bistatic range, and c stands for the speed of light.
Since R st and R sr are much larger than R tr , it can be assumed that the signal transmitted from the satellite is a plane wave rather than a spherical wave. Besides this assumption, a coefficient k is introduced to describe the relationship between R tr and R delta and can be expressed as When the bistatic angle is denoted as α and the angle between the receiver-to-target and receiver-to-satellite directions is represented as β, the coefficient k can be calculated by The geometric configuration of our radar used to detect the moving target (depicted as a ship) is shown in Fig. 1.
If the motion of the GEO satellite can be ignored, the locations of both receiver and illuminating satellite are fixed, and the delay difference between the direct and reflected signal is only related to the motion of the target. Based on the above approximations, after a cross-correlation between the received RC signal given by (5) and the reference signal reconstructed with the direct signal given by (4), the range-compressed signal can be written as where A is the signal amplitude and ρ is the envelope in the range-time dimension. The motion of the target can be usually expressed as a secondorder polynomial in a short period [34], i.e., a few seconds. In other words, the target moves translationally with a constant acceleration, and no high-order movement parameters such as jolt or rotational angular velocity need to be considered. Therefore, in terms of the radial acceleration a, the radial velocity v, and the initial target distance R 0 , R tr can be approximately expressed by After range compression, the one-dimensional rangecompressed signal is further divided and reformatted into a 2-D data array with a pulse repetition interval of 1 ms, which equals the period of the ranging code of the B3I component. To express the signal more clearly in a 2-D form, two parameters τ and t m are introduced here, which represent the fast-time and slow-time variables, respectively. The original time variable is the sum of these two parameters, i.e., t = t m + τ . The rearranged data in the 2-D form can be written as It is worth mentioning here that the envelope of the received signal in the azimuth-time dimension, which is related to the radiation pattern of receiving antenna, is neglected in our study. Furthermore, under the assumptions of the plane wave illumination and the second-order motion of the target, the expression of the IF signal after range compression can be found in (11), which provides a basic formula for subsequent studies in this article.

III. TARGET DETECTION AND PARAMETER ESTIMATION
For the received signal after range compression, the echo energy spread over different cells in the data array, as shown in (11). An energy-focusing process used to increase the SNR within one CPI is necessary. RCM and DCM corrections are the first two steps in the focusing procedure, and then the target detection and parameter estimation can be performed. The flowchart of our signal processing algorithm is shown in Fig. 2.
Before the description of the proposed algorithm, it is necessary to discuss the choice of CPI because this is one of the most important parameters in the subsequent signal processing and is also the main difference between our study and the conventional methods introduced in Section I. If the motion parameters of the target can be acquired before the detection, a longer CPI can be adopted to obtain a higher SNR because precise migration corrections can be achieved. Otherwise, the whole process of the target motion should be divided into several short frames to ensure that the target's movement over each frame can be expressed by low-order motion parameters. In contrast, a smaller CPI may reduce the obtained SNR. It can be found that there is an optimum value of CPI to acquire the largest SNR. In the traditional method, this value is found by progressively increasing the CPI until the SNR starts to decrease [12]. The purpose of this method is to improve the SNR as high as possible but at the expense of higher computational  complexity and less information about the target's motion. In our opinion, if the SNR is high enough to detect the target, it is also important to extract more information about the target's motion. Therefore, we decreased the CPI from the widely used 2.5 s until the echo signal was submerged by the noise. Using this strategy, we found that 1 s was the optimum value for the CPI, with which more information can be extracted to describe the target's movement process while the detection ability of the target does not deteriorate.
To better demonstrate the performance of our algorithm, the simulation results are also provided for each step of this process. The simulation parameters are given in Table I. It should be noted that the radial acceleration is set to be larger than the actual value, and hence, the echo signal in the azimuth-time domain can be more accurately expressed by the LFM signal.

A. RCM Correction With KT
Considering the envelope term ρ in (11), there is range walk (linear RCM), which is related to the radial velocity term v and range curvature (quadratic RCM) produced by the acceleration a; this is the first migration we need to correct. The RCM of the target in simulation is shown in Fig. 3(a) as range-compressed data.
Before implementing RCMC, we need to analyze the order of magnitude of the range curvature. The signal source used in our study is the B3I band of the BeiDou signal, and its PRN code has a chipping rate of 10.23 MHz. The range resolution achieved with this signal is c/2CR, i.e., 14.6 m, where CR is the signal chipping rate. Moreover, the CPI is always 1 s, and the acceleration of the ferry that we intend to detect is always less than 1 m/s 2 according to its actual movement parameter. Therefore, during one CPI, the range curvature is smaller than the size of one range-resolution cell if the movement can be represented by (10). Certainly, if the signal bandwidth becomes wider and the acceleration is high enough, the range curvature can exceed the size of one range-resolution cell. Nevertheless, it can be neglected in our study due to the low acceleration of the detected target.
Considering all factors mentioned above, the only task we need to handle is correcting the range walk produced by the first-order coefficient, i.e., the radial velocity. As mentioned before, we do not know the velocity information in advance. As a result, KT could deal with this problem efficiently. Data after range compression are first transformed into the range frequency domain, and the range walk in the envelope is converted into the phase term by performing DFT in the range time domain. When the range curvature is neglected, the obtained result can be written as follows: where f d = −kv/λ, γ a = −ka/λ, α = f r /f 0 , λ is the wavelength of the carrier signal, f r denotes the frequency in the range direction, and P (f r ) is the envelope in the range frequency domain after DFT.
Since the sampling rate in the azimuth-time domain or pulse repetition frequency (PRF) is fixed at 1 kHz, there might be a frequency alias of the Doppler frequency f d . Denoting the aliasing frequency as f d0 , we can obtain where n k is the aliasing coefficient. Then, the maximum non-aliasing radial velocity is λPRF/2, which equals 118.5 m/s. Thereupon, the n k is always 0 and f d equals f d0 because the radial velocity of the cargo ship is usually less than the maximum value. Therefore, we can solve the aliasing problem with minimal effort.
As shown in (12), both velocity and acceleration terms are coupled with the range frequency and cannot be separated. Here, we define a new slow-time variablet m to rescale the slow-time axis as a function of the range frequency in the following form: Substituting the original slow-time variable t m with the new termt m in (12), we obtain Because the range curvature is very small, as discussed before, we can neglect the term related to the acceleration given in (15). Moreover, for the B3I signal, because f r f 0 , so 1 + α ≈ 1, and hence, (15) can be rewritten as follows: After an inverse Fourier transform in the range-frequency domain, the received signal can be written as After RCM correction by KT, the envelope of data along the azimuth time domain is corrected to be flat. The simulation result is shown in Fig. 3(b).

B. DCM Correction With DA
After KT, there is still a quadratic component in the phase term, as given by the last term in (17). In other words, if we check the signal in the azimuth time domain, the received signal has a form of LFM signal and γ a is the chirp rate. To further increase the SNR of the received signal, we need to estimate this parameter and then eliminate it. Here, DA is employed to find the chirp rate. Note that, the coefficient k in the subsequent simulation part is set as 2 to more intuitively show the estimation of target distance and target radial motion parameters. When the receiver is located on the line connecting the target and the satellite, the above value of k can be satisfied. Also, in our simulation, the sampling rate of the IF signal is 62 MHz. Therefore, the size of one range bin is c/2f s with k = 2, i.e., about 2.4 m, where f s is the sampling rate of the IF signal. Before searching the chirp rate, we first need to select the index of the range cell that contains the focused energy of the target. As a result, we extract the column of azimuth time data, corresponding to τ = 2R 0 /c. The selected data can be written as follows: where e(t m ) denotes the noise. The frequency spectrum of this column of data is shown in Fig. 4(a), and we can clearly observe its LFM characteristics from the image. For the DA, we need to search a series of possible chirp rates γ a , and multiply the corresponding compensation phase term by the selected data. This phase term is defined by The compensated data will become a sinusoidal signal in the time domain or a single peak in the frequency domain if the estimated chirp rate is close enough to the actual chirp rate. For each obtained chirp rate, we transform the compensated signal into a frequency domain and calculate the power of the signal. The chirp rate pertaining to the largest energy is the expected result. This process can be written as follows: where argmax denotes an operation used to find the argument that gives the maximum value from a target function, and DFT represents the DFT operation.
To reduce the complexity of the search process, the algorithm can be divided into three steps. First, we search the chirp rate with a coarse step and achieve a rough estimation. Then, we perform a finer search with a smaller search step. Finally, we repeat the above procedures until we get an accurate estimation of the chirp rate. Specifically, the relationship between the modulation rate and radial acceleration is defined as γ a = −2a/λ, and the upper and lower limit of the radial acceleration searching range are set to be 2.5 and −2.5 m/s 2 , respectively, which are wide enough for a ferry target that we intended to detect. Therefore, the search range for the modulation rate varies from −21.4 to 21.4 Hz/s. Then, the initial search step Δγ is set to 1 Hz/s, which works for the first iteration. In each new iteration, the searching step is reduced by 10 times, and the modulation rate is searched around the previous peak location. In practical operation, there were five iterations in total, and after all the iterations were done, an accurate result can be obtained.
After the compensation with the correct chirp rate, the signal can be written as follows: After a DFT along the azimuth time, the final result can be written as follows: where f a denotes the azimuth frequency, sinc is the sampling function, and T expresses CPI. The simulation result of the selected data, as (22) shows, is plotted in Fig. 4(b). It is clear that the LFM-like signal is corrected to a single-frequency signal. The energy of the echo in the azimuth frequency domain has been focused, and the SNR of the peak increases from 8.7 to 16.5 dB. It should be noted that the actual clutter is more complicated than that adopted in our simulation, and the practical radial acceleration of the target can be much smaller than that used for simulation. As a result, the final SNR enhancement brought by DA will be lower than that obtained in this simulation.
It can be found that, based on the assumption that the range curvature can be neglected, the energy-focusing process is completed, the energy of the target is focused on the range cell corresponding to the initial location, and the frequency cell of the Doppler shift is related to the radial velocity. It is also worth noting that the initial location only represents the situation of this specific CPI.
Finally, the focused energy peak is compared with a preset threshold to detect the target. A peak in the RD domain, whose value is higher than the threshold, indicates a moving target. The simulation result of target detection in the RD domain is shown in Fig. 5. In this image, we can clearly see a peak that corresponds to the target, and its location is related to the target distance and Doppler shift of the reflected signal.

C. Radial Velocity Estimation
After RCM and DCM corrections, for those frames that contain the moving target, parameter estimation can be performed. Here, two methods were used. DA uses the Doppler shift information in the RD domain to estimate the radial velocity, and RTA uses the target distance information in the RD domain to estimate the radial acceleration.
Details of DA have already been introduced in Section III-A, so here we skip the detailed steps of this algorithm and focus on its processing results. DCM correction is completed by compensating for the second-order phase term with the estimated parameter γ a . After that, the peak location of echo in the azimuth frequency domain represents the Doppler shift f d0 . Then, the radial velocity can be estimated by Although DA is widely used for parameter estimation, it has several limitations. Because it is implemented on the data of each frame, the results only demonstrate the motion of the target for that specific period. It is difficult to assess the correctness of the output results of DA for one frame without comparison with others. Especially when the length of the target is relatively large in the range dimension, different parts of the target can produce different radial velocities. In this case, the RD map presents multiple Doppler peaks in one range resolution cell, and it is difficult to determine a specific peak to represent the kinematic characteristic of the target. In addition, during one specific CPI, the LFM characteristics can be ambiguous due to the short integration time, which means that the estimated acceleration is not reliable enough. Hence, for those RD maps with multiple peaks, we first list all of them and then choose the most reasonable one by comparing results with those acquired from adjacent maps. When it comes to the acceleration data, we simply discard them and keep the estimated results of radial velocity.

D. Radial Acceleration Estimation
Radial velocity is the derivative of the target distance, and radial acceleration is the derivative of radial velocity. Based on this idea, RTA uses target distance historic data to estimate the radial acceleration of the moving target. Specifically, radon transform (RT) is an integral transform and gives us a linear integral of one function along a line having a specific slope and a specific distance from the origin. As a result, RT can find all obvious slopes from this function. If we organize the target distance data as a polyline over time, the slopes will equal the radial velocities. Similarly, the slopes of a polyline of radial velocity equal the radial accelerations.
When we perform RT for the differentiated value of historic target distance data R tr , the estimated accelerations can be expressed as follows: where δ is the Dirac delta function and β is the rotation angle of the line used for integral. Furthermore, the distance between the line for the integral and the origin is set to 0 m. The radial acceleration a can be estimated with the angle β RTA is a process designed for the entire recorded data and only outputs obvious movement characteristics. At the same time, all the unreasonable data are ignored. The downside of RTA is that some movement parameters that last for a short period will be submerged by other long-term ones due to the principle of the integral. For the output results, it is also difficult to distinguish which phase of motion is related to the peaks. In addition, for some movement processes having many different phases, it is more reasonable to perform RTA on each phase separately.

IV. EXPERIMENTAL VERIFICATION
Our experiment was carried out on the riverbank of the Yangtze River in Wuhan. Fig. 6(a) shows the data acquisition geometry during the experiment, and Fig. 6(b) shows the azimuths and elevations of the BeiDou GEO satellites and the target by using the receiver's position as a reference. The receiver had one LHCP antenna pointing to the target area to acquire the weak echo signal and one RHCP antenna used for recording the direct signal. The experimental site is shown in Fig. 6(c), and please note that the antenna of the receiver was set with a depression angle of about 45°, which is the same as the elevations of GEO satellite C04. The azimuth of GEO C04 to the receiver is about 15°, and in the meanwhile, the bistatic angle is near this azimuth due to the distance between GEO C04 and the target being much larger than the target distance. Thus, according to (8), coefficient k is close to 2. Based on the above approximate conditions, we chose the BeiDou C04 satellite as the signal source for detection and took 2 as the value of k in the subsequent processing of the GEO C04 satellite data in this experiment. The moving target was a ferry going back and forth on the river and Fig. 6(d) shows a photograph of this ferry. The dimensions of this ferry were nearly 40 m by 10 m. During the measurement, a dataset of 100 s was recorded and further divided into 100 frames. Therefore, the period of each frame was 1 s. After RCM and DCM corrections, various peaks are more distinguishable in the RD domain, and hence all frames were processed into this form. The sounding result obtained from the second frame is shown in Fig. 7(a) as an example, where the unit of the horizontal axis is the target distance in meters and the unit of the vertical axis is converted into radial velocity in meters per second.
From Fig. 7(a), three peaks can be found. Since the reference signal is generated by the direct signal recorded by the HC channel, the peak located at the 0th range cell is exactly the direct signal component. Although the RC channel connects an LHCP antenna that can suppress the RHCP signal transmitted by the satellite directly, some of the energy of the direct signal can still leak into this channel, resulting in the range peak numbered as 0 in the image. On the right side of the peak generated by the direct signal, there is the echo energy scattered by the dock that is located near the receiver, as shown in Fig. 6(b). Because this ferry was stationary during the measurement, the peak generated by this target is located at the 0th frequency cell. The stationary target can be easily distinguished from the moving one with this characteristic in the frequency domain.
The peak on top was generated by the moving target that we intend to detect. Taking this scenario as an example, we could easily identify whether a moving target exists by selecting a reasonable peak from the recorded data in the RD domain. Moreover, the CLEAN method, in which the target with the strongest echo amplitude is removed in each iteration and finally the one with the weakest echo amplitude is remained [43], can be used to remove the energy produced by the stationary target and can further reduce the false alarm rate of this radar system. As shown in Fig. 7(a), the target is located within an area bounded by a velocity ranging from −5 to 5 m/s and a distance ranging from 50 to 300 m. Therefore, we utilize this pair of parameters to designate the boundaries of data to reduce the complexity of the fast Fourier transform calculation and RD processing because there is no need to process a frame in which many cells do not contain any echo energy. All 100 frames were analyzed, and the first 50 of them contain the target of interest. Through summing up these RD maps, we can get the target's trajectory in the RD domain, as shown in Fig. 7(b).
First, we qualitatively analyzed the entire process of the target's motion. When the data were recorded, the target ferry left the port and drove toward the left side of the river. Then, it turned around and drove toward the right side. The overall trend of radial velocity was to first accelerate and then decelerate. The actual velocity of the target was almost constant, but the radial velocity varied with the angle between the target's route and the observing direction. When the target turned around and the course was away from the receiver, the radial velocity reached its maximum. It is worth noting that, because the ferry drove away from the receiver, the value of radial velocity was always negative. For the target distance, it keeps increasing. In Fig. 7(b), it can be seen that both the radial velocity and the target distance of the target varied over time. Furthermore, the target's RD trajectory is displayed as a V-shaped curve, which is consistent with the actual movement of the ferry during the measurement.
Next, we conducted a quantitative analysis of this process by estimating radial velocity with DA and acceleration with RTA, respectively. Data from AIS are widely used as a reference for parameter estimation, but the CPI in our study was always 1 s, so AIS cannot provide a sufficiently reliable reference because of its unstable update rate. In this experiment, we collected real historic location data of the ferry by recording the position data from a GPS receiver on it, shaping a trajectory as Fig. 6(a) shows. The position data were recorded once per second; thus, this reference can provide a continuous comparison with the radar measurement. The observation direction of our radar system was recorded as well, so the target distance reference data can be computed by projecting the actual range of the target onto the observation direction. Then, radial velocity reference data can be obtained by differentiating the target distance data. As stated in the previous section, DA can output estimated results for each frame. Fig. 8(a) shows the radial velocity estimation results for the first 50 frames in the orange polyline. In the same subfigure, the reference data collected from a GPS receiver are shown in the blue polyline. It can be observed that the variation of the curve measured by our radar system is consistent with the actual movement curve of the target, although the error compared to the real data is neither stable nor negligible. That is because the GPS data have their own instability without optimization, and we do not further process the GPS data to improve their smoothness in this article.
For the target distance, the measured data of our system can be derived simply by converting the code phase of the peaks in the RD domain into meters. In Fig. 8(b), the target distance data estimated by our radar system and the reference data are plotted by an orange polyline and a blue polyline, respectively. It can be found that our estimated target distance data are accurate in comparison with the real data.
On the other hand, we implemented RTA on the derivative of target distance data to estimate the radial acceleration. The Radon map is shown in Fig. 9. RT outputs the angle of the rotated line for the integral, and the two most obvious angles are 85°and 98°, which can be converted to radial acceleration as −0.088 and 0.141 m/s 2 . Two estimated accelerations demonstrate the variation of radial velocity during the deceleration and acceleration stages in the whole movement process. Because GPS data cannot provide reference acceleration data, and the velocity data calculated from the target distance given by the GPS receiver vary dramatically, as shown in Fig. 8(a). Therefore, the GPS data cannot be used to calculate radial acceleration in this article. In such a case, the results of RTA were verified with the historic velocity data, which come from DA. During the first 30 s, the velocity decreased from −1.177 to −4.612 m/s, which can be represented by an acceleration of −0.106 m/s 2 . For the last 20 s, the velocity increased from −4.612 to −2.247 m/s, which can be represented by an acceleration of 0.118 m/s 2 . Due to the lack of real acceleration data, the results are of relatively low reliability, but at least produce a value that is useful as a reference in terms of the trend and magnitude of the radial velocity changes. From another point of view, as introduced in Section III, DA uses frequency domain data to estimate the velocity and RTA uses time domain data to estimate the acceleration. The results of these two algorithms are consistent with the actual condition and can also corroborate each other. Thus, the correctness and feasibility of both DA and RTA have been verified.
In addition, we processed the recorded signal with our proposed algorithms but with 3 s CPI to provide a comparison for further proving the practical advantage of our strategy. The RD trajectory generated by 3 s CPI processing is shown in Fig. 10. By comparing the two trajectories in Figs. 7(b) and 10, it is clear that the longer CPI does not provide us more information, and the generated results are of nearly the same SNR with our 1 s CPI algorithm. The reason for this SNR deterioration is that the motion in each segment is not well represented by a second-order polynomial. In other words, KT and DA cannot compensate for RCM and DCM as accurately as they can under the condition of 3 s CPI. Therefore, in order to estimate agile motion parameters, it is necessary to employ a small CPI for signal processing because we can interpret more information from the result with little loss of SNR, and a smaller CPI also brings less computational complexity, which can boost the processing speed.

V. MULTISATELLITE DATA FUSION
The most prominent merit of GNSS-based bistatic radar is the spatial diversity of the signal sources: there are hundreds of GNSS satellites in space. If an appropriate multi-satellite data fusion method can be developed, more information about the detected target could be extracted for target kinematic recovery and target identification. To achieve this objective, we proposed an algorithm to fuse echoes from multiple satellites and show the fusion result in one RD image.

A. Signal Source Selection
Our study focuses on the B3I component transmitted by BeiDou satellites. Specifically, the GEO satellites were selected among the available ones, and there were seven GEO satellites in operation when our experiment was carried out. The subsatellite point locations of BeiDou GEO satellites during our experiment are shown in Fig. 11, with the location of the experimental site.
According to the ICD, we know that GEO satellites transmit B3I signals. Moreover, we can see from Figs. 6(b) and 11 that the BeiDou C01, C04, and C59 satellites can satisfy our configuration better to maximize signal strength, as introduced in Section II. Thus, in this experiment, the echoes produced by the BeiDou C01, C04, and C59 satellites were chosen for data fusion.

B. Data Fusion
In practice, two methods can be adopted to fuse these data. One is coherent integration and the other is noncoherent integration. Coherent integration requires calculating and estimating the trajectories of two selected satellites in 3-D space and then compensating for the difference between their signal delays and Doppler shift to align their phase before plotting an RD map. Generally speaking, coherent integration requires maintaining the coherency of the signal phase when signals are integrated. In our experiment, various GEO satellites have different positions, resulting in different bistatic configurations when each of them was separately adopted as the transmitter. Since different bistatic configurations can produce different Doppler and phase shifts, the coherent integration of echo signals produced by multiple GEO satellites cannot be achieved without refined compensation for these echo signals. Therefore, we utilize the noncoherent integration method to integrate target returns to simplify the whole process, i.e., correcting the motion between RD trajectories and overlapping them. Since the SNR after the aforementioned signal processing is sufficiently high to locate the echo's peaks by visualization, we can extract the specific parameters of the echo's peaks from the BeiDou C01, C04, and C59 satellites in the same frame and compensate for the differences in the three echo peaks in both the range time domain and the azimuth frequency domain to achieve data fusion.
To fuse data, we first need to choose reference data. In practice, we selected the trajectory of the target when the BeiDou C04 satellite is used as the transmitter because this satellite satisfies our configuration best according to the previously introduced metrics. Theoretically, these trajectories could be simply fused by overlapping each other. However, due to the discrepancy between the actual bistatic configuration and our hypothesis, as mentioned in Section II about the geometric configuration, these trajectories have different variation trends. Therefore, considering the differences between the locations of the echo's peaks in RD maps, we need to compensate for these echoes before noncoherent integration in both the range domain and the Doppler domain. Below, we introduce our algorithm by presenting a procedure for the fusion of data when satellites C01 and C04 are separately used as the transmitter; this algorithm can also work for the fusion of data from other satellites.
For frame 1 in the data acquired by using BeiDou C01 satellite as the transmitter, the signal in the RD domain is represented as S RD_1_1 where subscript 1_1 denotes frame 1 in the data acquired by using BeiDou C01 satellite as the transmitter. For frame 1 in the data acquired by using the BeiDou C04 satellite as the transmitter, the signal is denoted as s 1_4 : where the subscript 1_4 represents frame 1 in the data acquired by using the BeiDou C04 satellite as the transmitter.
The difference between the two range cells is denoted as When the echo data produced by using BeiDou C04 satellite as the transmitter is selected as the reference, range cell compensation can be achieved by convoluting with a unit impulse signal at the location of rcd 1_4_1 in the time domain or by multiplying them in the frequency domain.
The difference between the two Doppler cells is denoted as The Doppler cell compensation can be achieved by multiplying the corresponding phase term in the time domain data or convoluting with a unit impulse signal at the location of dcd 1_4_1 in the frequency domain. After compensation in both the time and frequency domains, the two RD maps are fused, and the fused RD map can be represented as S RD_1 where A 1 = A 1_1 + A 1_4 is the new envelope after combination. This fusion method was applied to all 50 frames in practical processing for both C01 and C59, using data acquired with satellite C04 as the reference. The fused RD map is shown in Fig. 12.
Compared with Fig. 7(b), it can be found that the target's trajectory in the RD domain is clearer in both the latter half and the beginning part of the movement, and the radial velocity information is more visible. Noncoherent integration does improve the readability and continuity, but has limited performance against SNR enhancement, especially considering that there are only three satellites used for the data fusion. Note that, the performance of the noncoherent integration depends on the SNR of the original data. If the SNR of the original data is not sufficient, noncoherently integrating multiple data will worsen the fusion result. Although SNR enhancement in our study is limited by both available signal sources and the noncoherent integration algorithm, the data fusion still improves the quality of the sounding results in the RD map.

VI. CONCLUSION
This article investigates the possibility of detecting a moving target and estimating its parameters using BeiDou GEO satellites-based passive radar with a short integration time, i.e., 1 s. This purpose was achieved by exploiting KT and DA to correct RCM and DCM, respectively. This technique has advantages over other methods in terms of the integration time, which is of great importance, especially for real-time target detection with a GNSS-based passive radar system. Furthermore, we also introduced RTA, a novel algorithm used to estimate the radial acceleration of the target. Besides target detection and parameter estimation, a noncoherent integration method for data fusion with multiple satellites is also discussed, which provides more information in the RD domain. In addition, we recorded the location data of the target as the reference, and the correctness of our proposed signal processing method was verified. Target detection and parameter estimation are the prerequisites of further radar signal processing such as target tracking and target imaging, so we believe the technique introduced in this article is of great practical value.
However, further study is required on some topics. One is the localization of the target. All the parameters estimated with our algorithm are radial parameters, such as radial velocity and acceleration. In fact, the absolute location data, and the actual kinematics parameters of the target in motion could be obtained by solving the location equations expressed by coefficients from the parameters of at least three satellites. This is more like a cross-positioning technique and is theoretically achievable as long as the movement of each satellite can be precisely compensated. As a matter of fact, the aforementioned kinematics parameters recovery problem has been discussed for maritime targets in some recently published works [44], [45]. Another area for future study is the coherent integration of multiple illuminators. This is significant for a GNSS-based passive radar system because its performance needs to be continuously improved, while the power budget of a single illuminator and the SNR improvement achieved by noncoherent integration is limited.
Our future work will be in three areas. The first is to study the bistatic configuration in which non-GEO satellites are used as transmitters of opportunity, and an aircraft will be used as the detected target, which requires a finer and more complicated analysis of the dynamic delay difference, along with a more universal signal processing methodology. The second will focus on the coherent integration of the data from both GEO and non-GEO satellites. A final task is to collect echo data when non-GNSS satellites are used as transmitters of opportunity like communication satellites and broadcast satellites. The possibilities of the latter two are based on the prominent advantages of satellite-based passive radar systems, which are high availability and flexibility.