Moving Target Imaging Using GNSS-Based Passive Bistatic Synthetic Aperture Radar

: Current studies of global navigation satellite systems (GNSS)-based bistatic synthetic aperture radar (GNSS-SAR) is focused on static objects on land. However, moving target imaging is also very signiﬁcant for modern SAR systems. Imaging a moving target has two main problems. One is the unknown range cell migration; the other is the motion parameter estimation, such as the target’s velocity. This paper proposes a moving target imaging formation algorithm for GNSS-SAR. First, an approximate bistatic range history is derived to describe the phase variation of the target signal along the azimuth time. Then, a keystone transform is employed to correct the range cell migration. To address the motion parameter estimation, a chirp rate estimation method based on short-time Fourier transform and random sample consensus is proposed with high processing e ﬃ ciency and robust estimation errors in low signal-to-noise ratio scenes. The estimated chirp rate can calculate the target’s velocity. Finally, azimuth compression derivation is performed to accomplish GNSS-SAR imaging. A maritime experimental campaign is conducted to validate the e ﬀ ectiveness of the proposed algorithm. The two cargo ships in the SAR images have good accordance with the ground truth in terms of the target-to-receiver vertical distances along the range and the ships’ length along the cross-range.


Introduction
In the last decade, passive bistatic synthetic aperture radar (SAR) has received substantial attention from the research community. Because the transmitters are illuminators of opportunity, such as digital video broadcasting-terrestrial systems [1], frequency-modulated radio [2], GSM [3], and global navigation satellite systems (GNSS), only the receiver needs to be developed. Hence, a passive bistatic SAR has advantages of low cost, license-free, covertness, and no electromagnetic pollution. Among these opportunistic illuminators, GNSS-based bistatic SAR (GNSS-SAR) is a very promising research field. Besides the above advantages, GNSS-SAR systems have good system stability, short revisit time, and flexible configuration thanks to its global GNSS constellations, e.g., the global positioning system (GPS), global navigation satellite system (GLONASS), Galileo system, and the Beidou system.
The concept of GNSS-SAR was likely proposed for the first time in [4]. Then, the theoretical and experimental investigations were reported in [5,6], which confirmed the feasibility of the GNSS-SAR technique. To provide imagery, the derivation of an image formation algorithm has become a major challenge because GNSS-SAR has a flexible topology of bistatic acquisition. The transmitter is a GNSS satellite, while the receiver can be deployed on the ground, a moving vehicle, or an aircraft, resulting in various bistatic acquisition geometries. Several contributions to image formation algorithms for special bistatic SAR topologies can be found in recent papers, such as the range-Doppler algorithm [7,8], back-projection algorithm [9], range migration algorithm [10], and some improved algorithms [11][12][13]. The validation of the GNSS-SAR image formation algorithms has paved the way for GNSS-SAR application. For example, temporal land change monitoring based on GNSS-SAR was investigated in [14,15]. A GNSS-SAR imaging for the urban area with 360 • coverage was obtained in [16]. Moreover, considering the rough range resolution of GNSS-SAR, some papers focused on range resolution improvement [17][18][19][20].
The above studies on GNSS-SAR have been dedicated to static objects on land. However, to the best of the authors' knowledge, few contributions to moving target imaging have been presented in the GNSS-SAR field. With growing requirements for both civilian and military applications, moving target imaging has been a long-standing subject for modern SAR systems. Many investigations on moving target imaging have been conducted in the active SAR community [21][22][23][24]. In recent years, passive radar systems have been used for inverse SAR (ISAR) [25][26][27][28]. ISAR processing is a post-processing operation applied to the final SAR image [29]. To extend the range of GNSS-SAR applications, this paper concentrates on moving target imaging. The conventional GNSS-SAR image formation algorithms have two main steps, including range cell migration (RCM) correction and azimuth compression. For the stationary target, the RCM and azimuthal Doppler shifts of the echoes are dependent on the motion of the SAR platform (the receiver or the satellite). That is to say, RCM correction and azimuth compression can be achieved by the known motion of the platform. However, for the moving target, none of the above GNSS-SAR image formation algorithms can enable RCM correction and azimuth compression due to the unknown motion of the non-cooperative target. As a result, the moving target will be defocused in the SAR image. Therefore, moving target imaging has at least two main problems. One is the unknown RCM, and the second is the motion parameter estimation before azimuth compression.
In this paper, to address the main problems of moving target imaging, we select a relatively simple bistatic acquisition geometry shown in Figure 1. The receiver is fixed on the shore, and the GNSS satellite is assumed as stationary during the observation time (e.g., less than 120 s). A maritime moving target whose trajectory perpendicular to the line of sight (LOS) of the receiver antenna is selected as the object of interest so that a synthetic aperture can be formed. Based on the bistatic acquisition geometry, we first derive an approximate bistatic range history that can describe the phase variation of the target signal. To correct the RCM, a keystone transform [30] is employed without prior knowledge of the target's motion. Then, to deal with the parameter estimation problem, a chirp rate estimation method based on short-time Fourier transform (STFT) and random sample consensus (RANSAC) [31] is proposed with high processing efficiency and robust estimation errors. The estimated chirp rate can calculate the target's velocity. Finally, a derivation of azimuth compression is conducted to accomplish the moving target imaging. Two groups of real experimental data are obtained to validate the effectiveness of the proposed method, and then the experimental results are reported and discussed in detail.
The paper is organized as follows. An approximate bistatic range history is derived in Section 2. Then, Section 3 describes the image formation algorithm. Section 4 provides the experimental setup and results. A further discussion about the experimental results is presented in Section 5, and Section 6 concludes this paper. The paper is organized as follows. An approximate bistatic range history is derived in Section 2. Then, Section 3 describes the image formation algorithm. Section 4 provides the experimental setup and results. A further discussion about the experimental results is presented in Section 5, and Section 6 concludes this paper.

Approximate Bistatic Range History
In this section, the bistatic acquisition geometry of maritime moving target imaging is introduced to derive an approximate bistatic range history that can describe the phase variation of the target signal for the next image formation algorithm. At the same time, the maximum error between the approximate bistatic range and the real bistatic range will be analyzed. Figure 1 shows the bistatic acquisition geometry. The receiver on the shore has a pair of channels: the reference channel (RC) steers the antenna toward the sky to collect the direct GNSS signals as reference, and the surveillance channel (SC) points the antenna toward the sea surface to record the reflected GNSS signals from the moving target.
For the sake of simplicity, the receiver and the target are located at the same height. The simple side view of the bistatic acquisition geometry is illustrated in Figure 2, where the satellite (S), receiver (R), and target (T) constitute a plane. It is assumed that the satellite in space is approximated as stationary during the observation time (e.g., less than 120 s). The bistatic range history of the moving target during the observation time is defined as: where ∈ − , is the azimuth time (or slow time), T is the entire observation time, is the transmitter-to-target range, is the target-to-receiver range, and is the baseline range between the satellite and the receiver. Note that the bistatic range history in (1) is accurate without any approximation. Nevertheless, it is not convenient for the next derivation of the moving target image formation algorithm. Hence, an approximate bistatic range history is introduced here.

Approximate Bistatic Range History
In this section, the bistatic acquisition geometry of maritime moving target imaging is introduced to derive an approximate bistatic range history that can describe the phase variation of the target signal for the next image formation algorithm. At the same time, the maximum error between the approximate bistatic range and the real bistatic range will be analyzed. Figure 1 shows the bistatic acquisition geometry. The receiver on the shore has a pair of channels: the reference channel (RC) steers the antenna toward the sky to collect the direct GNSS signals as reference, and the surveillance channel (SC) points the antenna toward the sea surface to record the reflected GNSS signals from the moving target.
For the sake of simplicity, the receiver and the target are located at the same height. The simple side view of the bistatic acquisition geometry is illustrated in Figure 2, where the satellite (S), receiver (R), and target (T) constitute a plane. It is assumed that the satellite in space is approximated as stationary during the observation time (e.g., less than 120 s). The bistatic range history of the moving target during the observation time is defined as: where t m ∈ − T 2 , T 2 is the azimuth time (or slow time), T is the entire observation time, R t is the transmitter-to-target range, R r is the target-to-receiver range, and R b is the baseline range between the satellite and the receiver. Note that the bistatic range history in (1) is accurate without any approximation. Nevertheless, it is not convenient for the next derivation of the moving target image formation algorithm. Hence, an approximate bistatic range history is introduced here. Because the transmitter-to-target range is in an order of tens of thousands of kilometers, much larger than the target-to-receiver range, line SR and line ST in Figure 2 can be regarded as parallel. The bistatic range history defined in (1) can be approximated as: Because the transmitter-to-target range is in an order of tens of thousands of kilometers, much larger than the target-to-receiver range, line SR and line ST in Figure 2 can be regarded as parallel. The bistatic range history defined in (1) can be approximated as: where θ is the bistatic angle. However, the bistatic angle cannot be measured directly and varies with the target's motion. The three-dimensional view of the bistatic acquisition geometry is presented in Figure 3a for the bistatic angle analysis. By adding some auxiliary lines (dashed lines) in the solid geometry, we have: Remote Sens. 2020, 12, x 5 of 21 Moreover, if the target in Figure 3b moves towards the opposite direction, likewise, the bistatic range is: In summary, before the target crosses the LOS of the receiver antenna, if the target and the nadir point of the satellite are on the same side in Figure 3b with respect to the receiver, the bistatic range expression is (8); otherwise, it is (9). This peculiarity will be used in Section 5.2.  To evaluate the error between the approximate bistatic range and the real bistatic range, simulation calculation is conducted. The parameters in Table 1 are used to calculate the real bistatic range in (1) and the approximate bistatic range in (8) or (9). Note that the local satellite azimuth angle in Table 1 is equal to the satellite azimuth angle for convenience. Figure 4 shows the maximum errors between the approximate bistatic ranges and the real bistatic ranges versus the target-to-receiver vertical ranges for various satellite elevation angles and azimuth angles. We can observe that the maximum errors become large gradually with the increase of satellite elevation angles, azimuth angles, and target-to-receiver vertical ranges. The maximum errors do not exceed 1 m within the range of 5000 m, which is acceptable. Therefore, the approximate bistatic range history can be safely used for the next moving target image formation algorithm.

Parameters Values
Coordinates of the receiver (0, 0, 0) Satellite-to-receiver range (Rb) 20,000 km Satellite elevation angle (α) 0°, 30°, 45°, 60° Local azimuth angle (az) 0°, 30°, 45°, 60° Radar antenna beam angle (θbw) 10° Target-to-receiver vertical range (Rs) 500-5000 m With algebraic simplification, the bistatic angle is derived as: where α is the satellite elevation angle and β is a vertex angle constituted by the receiver, target, and nadir point of the satellite. The top view of the bistatic acquisition geometry in Figure 3b is presented to analyze β. In Figure 3b, a point target T is moving from the initial position T0 to the right side. The trajectory of the point target as the synthetic aperture is perpendicular to the LOS of the receiver antenna (assuming that the target moves within the antenna footprint during the whole observation time). Some angles can be obtained as: and where β and ε are the vertex angles, ϕ is the squint angle between R r and R s , R s is the target-to-receiver vertical range (i.e., the closest range), and az is named as local satellite azimuth angle that is the difference between the satellite azimuth angle and the azimuth angle of the LOS of the receiver antenna Remote Sens. 2020, 12, 3356 5 of 20 (Section 4 will give the computing equation according to the field test). Nevertheless, ε is so small that β ≈ β. Substituting (6) into (4), we have the bistatic angle as: where t m ∈ 0, 2L v , v is the target's velocity, and L = |x r − x 0 | is the half length of the synthetic aperture. Equation (7) shows that the variation of the bistatic angle θ is related to the target's motion. Substituting (7) into (2), we obtain the approximate bistatic range history as: Moreover, if the target in Figure 3b moves towards the opposite direction, likewise, the bistatic range is: In summary, before the target crosses the LOS of the receiver antenna, if the target and the nadir point of the satellite are on the same side in Figure 3b with respect to the receiver, the bistatic range expression is (8); otherwise, it is (9). This peculiarity will be used in Section 5.2.
To evaluate the error between the approximate bistatic range and the real bistatic range, simulation calculation is conducted. The parameters in Table 1 are used to calculate the real bistatic range in (1) and the approximate bistatic range in (8) or (9). Note that the local satellite azimuth angle in Table 1 is equal to the satellite azimuth angle for convenience. Figure 4 shows the maximum errors between the approximate bistatic ranges and the real bistatic ranges versus the target-to-receiver vertical ranges for various satellite elevation angles and azimuth angles. We can observe that the maximum errors become large gradually with the increase of satellite elevation angles, azimuth angles, and target-to-receiver vertical ranges. The maximum errors do not exceed 1 m within the range of 5000 m, which is acceptable. Therefore, the approximate bistatic range history can be safely used for the next moving target image formation algorithm. Table 1. Parameters for the real bistatic range and the approximate bistatic range calculation.

Moving Target Image Formation Algorithm
The block diagram of the proposed moving target image formation algorithm for GNSS-SAR is presented in Figure 5. Each processing block will be described in the following.

Moving Target Image Formation Algorithm
The block diagram of the proposed moving target image formation algorithm for GNSS-SAR is presented in Figure 5. Each processing block will be described in the following.

Signal Synchronization and Signal Model
In the passive bistatic SAR system, a direct signal is usually employed as the reference signal for the target signal compression. However, the signal-to-noise ratio (SNR) of the direct signal output from the RC antenna can be as low as −30 dB [5]. For this reason, it is critical to track the parameters

Signal Synchronization and Signal Model
In the passive bistatic SAR system, a direct signal is usually employed as the reference signal for the target signal compression. However, the signal-to-noise ratio (SNR) of the direct signal output from the RC antenna can be as low as −30 dB [5]. For this reason, it is critical to track the parameters of the direct signal, including the code phase, carrier phase, Doppler, navigation bits, etc. These parameters can generate a noise-free replica of the direct signal. The above processing is called signal synchronization and is very significant for the GNSS-SAR system. To acquire and track the direct signal, a software receiver is exploited. In our previous work, the software receiver has been successfully applied for deformation monitoring based on GNSS-R technology in [32], which implies that the software receiver is reliable and stable.
With the tracked parameters, the local reference signal is generated for the next range compression for the target signal. Here, taking GPS L1 signal as an example and ignoring the constant phase and amplitude terms, we can express the local reference signal as radar data formatting: where t ∈ [0, PRI] is the fast time, PRI means the pulse repetition interval (i.e., the duration time of the C/A code), t m ∈ − T 2 , T 2 is the azimuth time, T is the entire observation time, C(•) is the C/A code whose duration time is 1 ms, D(•) is the navigation data, and e(•) denotes the complex carrier. τ d , f d and ϕ d (t m ) are the code phase delay, carrier frequency, and carrier phase term of the local reference signal, respectively. Their values are dependent on the relative position and movement between the satellite and receiver. Moreover, τ d and ϕ d (t m ) contain the total delay and phase errors, respectively, induced by the atmospheric factors (i.e., troposphere and ionosphere delay) and the receiver errors (i.e., clock cycle slip and local oscillator drift). Likewise, assuming a point target case, we can express the target signal in the SC as: where τ r , f r , and ϕ r (t m ) are the code phase delay, carrier frequency, and carrier phase term of the target signal, respectively.

Range Compression
To compress the target signal along the fast time and preserve the exponential phase term for the next azimuth-time processing, range compression operation is implemented. By cross-correlating the complex conjugate of (10) with (11) for every azimuth time, we have range compression as: where * denotes the complex conjugate. ∆ f (t m ) and ∆ϕ(t m ) are the instantaneous difference between the local reference signal and the target signal with respect to carrier Doppler frequency and phase, respectively. Note that: • The reference signal and the target signal have the same navigation data within the range of 6000 km [33]. • ∆ϕ(t m ) has eliminated the total phase errors, respectively, induced by the atmospheric factors and the receiver errors due to the similar atmospheric factors and the shared oscillator in the RC and SC [34].
can be neglected considering the low Doppler frequencies induced by the moving target during the PRI [34]. After integration, we can rewrite (12) as: where CF(•) represents the envelope of the cross-correlation function, R bi (t m ) is the bistatic range history, c is the speed of light, and f c is the central carrier frequency.

Range Cell Migration Correction
From (13), we can observe that the varying bistatic range gives rise to both RCM and Doppler shifts along the azimuth time. The Doppler shifts can be used for motion parameter estimation and imaging later. For the RCM, if the range resolution is coarse and the target is not far from the receiver, there is no need for an additional RCM correction in the algorithm, as shown in Figure 6a. For instance, the best range resolution of GPS-based SAR is about 150 m in a quasi-monostatic configuration [20]. Otherwise, the RCM should be tackled and compensated.
• has eliminated the total phase errors, respectively, induced by the atmospheric factors and the receiver errors due to the similar atmospheric factors and the shared oscillator in the RC and SC [34].
• can be neglected considering the low Doppler frequencies induced by the moving target during the PRI [34].
After integration, we can rewrite (12) as: where • represents the envelope of the cross-correlation function, is the bistatic range history, is the speed of light, and is the central carrier frequency.

Range Cell Migration Correction
From (13), we can observe that the varying bistatic range gives rise to both RCM and Doppler shifts along the azimuth time. The Doppler shifts can be used for motion parameter estimation and imaging later. For the RCM, if the range resolution is coarse and the target is not far from the receiver, there is no need for an additional RCM correction in the algorithm, as shown in Figure 6a. For instance, the best range resolution of GPS-based SAR is about 150 m in a quasi-monostatic configuration [20]. Otherwise, the RCM should be tackled and compensated.  Let us take the instant time when the moving target is crossing the LOS of the receiver antenna as the reference time (here, it is assumed to be t m = 0, corresponding to the target at the imaging scene center). Taking (8) for derivation and applying a second-order Taylor series expansion, we can obtain R bi (t m ) as: where t m ∈ − L v , L v . Equation (14) shows that the RCM is related to the target's velocity. However, in practice, the velocity of the non-cooperative target is unknown. Here, a keystone transform is used to correct the RCM without prior knowledge of the target's velocity. To apply keystone transform, the range-compressed target signal in (13) is put in the range frequency and azimuth time domain by Fourier transform (FT), which is written as: where f r is the range frequency component and CF( f r ) is the envelope of the target signal in the range frequency dimension. By inserting (14) into (15), we have the exponential phase term of the signal: (16), we can see that the coupling of the range frequency and the azimuth time in the first and the second phase term results in range curvature and range walk, respectively. The range walk can be removed if we rescale the azimuth time for each range frequency by keystone transform as: Substituting (17) into (16), the rescaled phase term in the new ( f r , τ m ) domain is: where λ = c f c is the central wavelength of the carrier. Because the GNSS signal is not designed for SAR application, the narrow band of GNSS ranging code makes f r f c , and we have (18) can be rewritten as: As shown in (19), the first square brackets correspond to the quadratic phase term of the target signal along the azimuth time, while the second square brackets indicate that the RCM is compensated into the constant bistatic range. The above procedure allows the coupling effects of the range walk to be removed, but the residual range curvature is still present, as shown in Figure 6b. However, the effect of range curvature is expected to be smaller than a one range resolution cell since Galileo E5a/b or GPS L5 signals with a bandwidth of 10.23 MHz can provide the best range resolution of 15 m in a quasi-monostatic configuration [6].
After the RCM correction, a range inverse Fourier transform (IFT) is conducted for (15) to derive the range-compressed target signal as: where R c = R s (1 + cosαcosaz) is the constant bistatic range. The cross-correlation envelopes of the target signal along the azimuth time are now compensated to the same range cell. Hence, the target-to-receiver vertical range R s can be estimated by: The target-to-receiver vertical range will be used later in Sections 3.4 and 3.5.

Motion Parameter Estimation
Because the target's motion induces the modulation of the target signal along the azimuth time, the motion parameter, such as the target's velocity, should be obtained before azimuth compression. Otherwise, the moving target will be smeared and shifted in a static SAR image. In this paper, we refer to a target moving at almost constant velocity, such as ships sailing at cruising speed. Therefore, the target velocity estimation issue will be tackled in this subsection.
The Doppler shifts of the target signal can be obtained from the first derivative of the exponential phase term of (20) versus the azimuth time as: where γ is the chirp rate and f 0 is the Doppler centroid. Equation (22) indicates that the target signal along the azimuth time is a linear frequency modulation (LFM) signal and the target's velocity is related to the constant chirp rate. The velocity estimation issue can be transformed into the chirp rate estimation of the LFM signal. Time-frequency analysis techniques are effective methods for chirp rate estimation, such as Radon-Wigner transform [35], chirp-Fourier transform [36], and fractional Fourier transform [37].  Figure 7 shows the flowchart of the processing chain.
Because the target's motion induces the modulation of the target signal along the azimuth time, the motion parameter, such as the target's velocity, should be obtained before azimuth compression. Otherwise, the moving target will be smeared and shifted in a static SAR image. In this paper, we refer to a target moving at almost constant velocity, such as ships sailing at cruising speed. Therefore, the target velocity estimation issue will be tackled in this subsection.
The Doppler shifts of the target signal can be obtained from the first derivative of the exponential phase term of (20) versus the azimuth time as: where is the chirp rate and is the Doppler centroid. Equation (22) indicates that the target signal along the azimuth time is a linear frequency modulation (LFM) signal and the target's velocity is related to the constant chirp rate. The velocity estimation issue can be transformed into the chirp rate estimation of the LFM signal. Time-frequency analysis techniques are effective methods for chirp rate estimation, such as Radon-Wigner transform [35], chirp-Fourier transform [36], and fractional Fourier transform [37]. These methods usually require a one-dimensional search and projection integration for the chirp rate estimation. Let be the number of time samples. The computational costs of these time-frequency methods are in the order of . However, long aperture time is essential in the GNSS-SAR to improve the azimuth resolution and the SNR, which leads to these time-frequency analysis methods to be very time-consuming. To counteract this problem, a chirp rate estimation method based on STFT and RANSAC is proposed. The computational cost of the proposed method is in the order of + because STFT is implemented on a computer using fast Fourier transform (FFT) with the order of and RANSAC needs times of iterations. Figure 7 shows the flowchart of the processing chain.
where * • is the complex conjugate of the analyzing window. Figure 8 shows two joint timefrequency distributions produced by the STFT, corresponding to two cargo ships in the field test. We can observe that two straight lines are displayed with width because the target signal in the real scene where win * (•) is the complex conjugate of the analyzing window. Figure 8 shows two joint time-frequency distributions produced by the STFT, corresponding to two cargo ships in the field test. We can observe that two straight lines are displayed with width because the target signal in the real scene is the sum of echoes from multiple scattering points on the target surface. Moreover, the slopes of the straight lines are, in fact, the chirp rates of two cargo ships.
To estimate the slope values, i.e., fitting lines, the straight lines should be extracted first. A suitable threshold T is used to separate the time-frequency sampling points of the straight line (target points) from the background noise sampling points (noise points) by: where A is the magnitude value of the sampling point in the joint time-frequency distribution. These sampling points construct a sampling set S = (t 1 , f 1 ), (t 2 , f 2 ), . . . , (t n , f n ) . Figure 9 shows two binary images generated by the sampling sets of two cargo ships.
Remote Sens. 2020, 12, x 11 of 21 is the sum of echoes from multiple scattering points on the target surface. Moreover, the slopes of the straight lines are, in fact, the chirp rates of two cargo ships. To estimate the slope values, i.e., fitting lines, the straight lines should be extracted first. A suitable threshold T is used to separate the time-frequency sampling points of the straight line (target points) from the background noise sampling points (noise points) by: where is the magnitude value of the sampling point in the joint time-frequency distribution. These sampling points construct a sampling set = , , , , … , , . Figure 9 shows two binary images generated by the sampling sets of two cargo ships.

RANSAC and Median Filter
Due to the low SNR, it is impossible to ensure without the noise points shown in Figure 9. The number of noise points may even be more than that of the target points. Therefore, traditional least square method (LSM) may fail to be applied in the slope value estimation (Section 4.2 shows an  To estimate the slope values, i.e., fitting lines, the straight lines should be extracted first. A suitable threshold T is used to separate the time-frequency sampling points of the straight line (target points) from the background noise sampling points (noise points) by: where is the magnitude value of the sampling point in the joint time-frequency distribution. These sampling points construct a sampling set = , , , , … , , . Figure 9 shows two binary images generated by the sampling sets of two cargo ships.

RANSAC and Median Filter
Due to the low SNR, it is impossible to ensure without the noise points shown in Figure 9. The number of noise points may even be more than that of the target points. Therefore, traditional least square method (LSM) may fail to be applied in the slope value estimation (Section 4.2 shows an

RANSAC and Median Filter
Due to the low SNR, it is impossible to ensure S without the noise points shown in Figure 9. The number of noise points may even be more than that of the target points. Therefore, traditional least square method (LSM) may fail to be applied in the slope value estimation (Section 4.2 shows an experimental example). RANSAC can perform parameter estimation in low SNR because it can address sampling sets involving many wrong data points by iteration. The general processing steps of RANSAC were stated in [31]. In our case, the model is (22), and RANSAC is employed to fit an optimal straight line from S and derive the estimated slope value. Moreover, several optimal slope values may be obtained because of the straight lines with width in Figure 9. For robust estimation, after RANSAC operation, a median filter is used to get the suitable slope value with a high confidence level. Inliers are denoted as the target points that can be expressed by the model, while outliers are the noise points. The steps of RANSAC iteration are demonstrated below: Step 1: Let K be the maximum iteration number and generate S by STFT.
Step 2: Select two points randomly from S and construct the line via (22). Note that the process of random selection has some limitations that can improve the robustness:

•
Because the chirp rate in (22) is always negative, two points with the positive slope value are not considered.

•
The maximum considered chirp rate value can be set to limit the point selection region.

•
Because of the long observation time, two points with respect to the time axis should be selected as far away as possible so that the line can fit more sampling points.
Step 3: For every point (t i , f i ), the distance of point (t i , f i ) to the constructed line in Step 2 is calculated by: where |•| represents the absolute value. If d i ≤ th where th is the tolerance threshold, point (t i , f i ) is regarded as an inlier.
Step 4: Count the number of inliers as N k . If N k ≥ N s where N s is the minimum number of inliers that agrees with the line, the slope value of the line is put into the consensus set.
Step 5: Let k = k + 1. If k ≤ K, go back to Step 2. Otherwise, proceed to Step 6.
Step 6: A median filter is used to select the suitable slope value from the consensus set. Finally, the target's velocity is calculated via the target-to-receiver vertical range in (21) and the chirp rate equation in (22).

Azimuth Compression
The principle of azimuth compression is to remove the modulation phase term of the target signal through the designed azimuth matched filter in the range and Doppler frequency domain. Therefore, a derivation of azimuth matched filter is given in this subsection.
To design the azimuth matched filter, we should first obtain the analytical solution of the target signal in the azimuth frequency domain. For simplicity, one scattering point is considered here. The range-compressed target signal in (20) is transformed into the range and azimuth frequency domain by FT, which is written as: where f d is the azimuth frequency component and R bi (τ m ) is expressed by (8) or (9). However, it is hard to directly obtain the analytical solution of (26). The principle of stationary phase [38] is employed to address (26). After some lengthy algebra, we have: or where σ represents the magnitude of the Doppler spectrum, which is not important. From the previous time-frequency analysis, we know that the scattering points on the target surface experience the same Doppler shifts. Hence, Doppler spectra of N scattering points can be extended as: ) or where (•) represents linear superposition of the N scattering points' echoes, L n = |x r − x n |, and x n is the azimuth position. We can see that there are three exponential phase terms in (29) and (30). The first term is a modulation phase term that is the same for every scattering point. The second term is a linear phase term that contains the azimuth position of the scattering point. The last term is a constant phase term that has no effect on imaging. Obviously, the azimuth matched filter is a complex conjugate of the modulation term in (29) or (30), derived as: Equations (31) and (32) show that which azimuth matched filter is chosen depends on the target's moving direction. For example, (31) corresponds to the target moving towards the right side in Figure 3b, while (32) corresponds to the opposite direction. In other words, combining the bistatic acquisition geometry and the azimuth matched filter, we can determine the target's moving direction. This peculiarity will be discussed later in Section 5.2.
The azimuth matched filter generation requires parameters including the satellite elevation angle, the satellite azimuth angle, the azimuth angle of the LOS of the receiver antenna, the target-to-receiver vertical range, and the target's velocity. The first three parameters can be obtained through ephemeris data and measurement. The target-to-receiver vertical range and the target's velocity have been obtained in Sections 3.3 and 3.4, respectively. Then, to remove the modulation phase term, (29) is multiplied by (31), while (30) is multiplied by (32). An azimuthal IFT is conducted to accomplish the azimuth compression. The whole azimuth compression is derived as: where B = 2v×sin θ λ is the Doppler bandwidth, θ is the maximum squint angle, and t n m = L n v is the instant slow time when the n th scattering point crosses the antenna point direction.

Experimental Setup
A maritime experimental campaign was carried out at 16:25 on 16 May 2019. A GPS L1 signal was chosen as the non-cooperative illuminator of opportunity. The receiver was fixed at Cyberport Waterfront Park in Hong Kong. Figure 10a shows the receiver front end that has a pair of channels to collect the direct and reflected signals. Two different types of antennas shown in Figure 10b connect with the receiver front end. The circular antenna was a commercial off-the-shelf right-hand circular polarization (RHCP) antenna for direct signal collection, while the square antenna was a custom left-hand circular polarization (LHCP) antenna with a 12 dB gain for reflected signal collection. One can obtain more details on the receiver front end and the software receiver in Section 3 of [32]. An automatic identification system (AIS) [39] was exploited to provide voyage-related information of the ships as the ground truth. Table 2 lists some essential parameters for the imaging processing. From the zenith, Figure 11 presents the top view of the real scene where the azimuth angles of both the satellite and the LOS of the receiver antenna (clockwise from north) are marked. Hence, the local satellite azimuth angle defined in (6) is calculated by: where Sat_Az denotes the satellite azimuth angle and LOS_Az is equal to the azimuth angle of the LOS of the receiver antenna with the subtraction of 180 • . We selected two cargo ships shown in Figure 12 as the main objects. Table 3 gives the AIS information of the ships and the selected GPS satellites including the pseudo random noise code (PRN) numbers, the elevation angles, and the azimuth angles.
To validate the effectiveness of the proposed image formation algorithm, the collected raw data were addressed in MATLAB (MathWorks, Natick, MA, USA) on a computer with Intel i7-7700 3.60 GHz processors and 64 GB RAM. The results will be described later in this work.
where _ denotes the satellite azimuth angle and _ is equal to the azimuth angle of the LOS of the receiver antenna with the subtraction of 180°. We selected two cargo ships shown in Figure  12 as the main objects. Table 3 gives the AIS information of the ships and the selected GPS satellites including the pseudo random noise code (PRN) numbers, the elevation angles, and the azimuth angles. To validate the effectiveness of the proposed image formation algorithm, the collected raw data were addressed in MATLAB (MathWorks, Natick, MA, USA) on a computer with Intel i7-7700 3.60 GHz processors and 64 GB RAM. The results will be described later in this work.
(a) (b)  The tolerance threshold (PRF/T) × 3 The proportion of the minimum number of inliers 7.5%      Figure 13 shows the fitting results of two cargo ships via RANSAC, LSM, and the ground truth from AIS. In Figure 13a, we can find that the number of noise points seem to be more than that of the target points. However, RANSAC can still provide a good fitting line (red line) compared with the ground truth (black line), which implies that the estimated velocity is close to the AIS velocity. However, LSM shows a fitting line (green dashed line) with a positive slope value. As described in (22), the chirp rate is always negative. Therefore, LSM has failed to be applied in the velocity estimation. In Figure 13b, there is a small angle between the red line and the black line. The main reason is that the fluctuations of the target complex reflectivity influence the distribution of the sampling points. Particularly, the echoes whose Doppler frequencies are near zero have the strongest scattering effect shown in Figure 8b, resulting in more sampling points than that from other Doppler  Figure 13 shows the fitting results of two cargo ships via RANSAC, LSM, and the ground truth from AIS. In Figure 13a, we can find that the number of noise points seem to be more than that of the target points. However, RANSAC can still provide a good fitting line (red line) compared with the ground truth (black line), which implies that the estimated velocity is close to the AIS velocity. However, LSM shows a fitting line (green dashed line) with a positive slope value. As described in (22), the chirp rate is always negative. Therefore, LSM has failed to be applied in the velocity estimation. In Figure 13b, there is a small angle between the red line and the black line. The main reason is that the fluctuations of the target complex reflectivity influence the distribution of the sampling points. Particularly, the echoes whose Doppler frequencies are near zero have the strongest scattering effect shown in Figure 8b, resulting in more sampling points than that from other Doppler frequency regions. The estimated chirp rate and velocity are listed in Table 4. Considering the fluctuations of the target complex reflectivity, the errors between the estimated velocities and the AIS velocities are reasonable.     Figure 14a,b, respectively. The bright spots corresponding to two cargo ships provide good characteristics in terms of signal to background ratio. The range values of the bright spots have good accordance with the vertical range in Table 3, at less than one range resolution cell of the GPS L1 signal. Due to the high azimuth resolution of GNSS-SAR, Figure 14a,b also demonstrate the ship length measurement results along the cross-range: WAN HAI 506 is about 254 m, and WAN HAI 313 is about 204 m. The measurement results are very close to the real targets' length in Table 3. Moreover, Figure 14b indicates two bright spots located at different ranges, corresponding to the cargo ship of interest (ellipse box) and the yacht (square box) in Figure 12b, which further confirms the effectiveness of the proposed image formation method. For comparison, Figure 14c

Discussion
In this section, the mean square error (MSE) of the proposed chirp rate estimation method and the target's moving direction determination will be discussed.

Mean Square Error of the Chirp Rate Estimation Method
To quantitatively evaluate the chirp rate estimation performance of the proposed motion parameter estimation method in Section 3.4, Monte Carlo trials were conducted. At the same time, LSM and Cramer-Rao lower bound (CRLB) [40] were employed for comparison. Table 5 lists the detailed parameters for the simulations. Gaussian white noise as the background disturbance was added to the target signals, while no sea clutter was considered. Note that the parameters of RANSAC are the same as that in Table 2. Figure 15 shows the curves of the MSEs of the chirp rates versus  Table 4, (b) is WAN HAI 313 with the estimated velocity in Table 4, (c) is WAN HAI 506 with target's velocity = 3 m/s, and (d) is WAN HAI 313 with target's velocity = 3 m/s.

Discussion
In this section, the mean square error (MSE) of the proposed chirp rate estimation method and the target's moving direction determination will be discussed.

Mean Square Error of the Chirp Rate Estimation Method
To quantitatively evaluate the chirp rate estimation performance of the proposed motion parameter estimation method in Section 3.4, Monte Carlo trials were conducted. At the same time, LSM and Cramer-Rao lower bound (CRLB) [40] were employed for comparison. Table 5 lists the detailed parameters for the simulations. Gaussian white noise as the background disturbance was added to the target signals, while no sea clutter was considered. Note that the parameters of RANSAC are the same as that in Table 2. Figure 15 shows the curves of the MSEs of the chirp rates versus different SNR levels, and in each case, 100 Monte Carlo simulations are done. In Figure 15, if the SNR is more than −50 dB, the MSEs of LSM are lower than that of RANSAC and even close to the CRLBs. However, the estimation performance of LSM will reduce sharply with the SNR decrease because a large number of outliers in the time-frequency sampling set lead to the fitting process being distorted. Compared with LSM, RANSAC demonstrates advantages for line-fitting in the presence of outliers. The MSEs of RANSAC in Figure 15 can keep stable and approach a relatively low level when the SNR is lower than −55 dB. In practical application, it is more inclined to have good estimation performance when the input SNR is low. Therefore, the proposed motion parameter estimation method is suitable for the GNSS-SAR application. is lower than −55 dB. In practical application, it is more inclined to have good estimation performance when the input SNR is low. Therefore, the proposed motion parameter estimation method is suitable for the GNSS-SAR application. Target-to-receiver vertical range 1000 m Figure 15. MSEs of chirp rates by using LSM and random sample consensus (RANSAC).

Target's Moving Direction Determination
As deduced in Section 3.5, combining the bistatic acquisition geometry and the azimuth matched filter, we can determine the target's moving direction. Take WAN HAI 506 for example. First, we compare the SAR images produced by two azimuth matched filters defined in (31) and (32). Figure  14a is produced by using (32), in which we can see that the target signal energy is concentrated together. Figure 16 indicates the SAR image with (31). Due to mismatch during the azimuth compression, the target signal energy in Figure 16 has been split into two parts. Therefore, the correct azimuth matched filter is (32). According to the derivation in Sections 2 and 3.5, (32) corresponds to

Target's Moving Direction Determination
As deduced in Section 3.5, combining the bistatic acquisition geometry and the azimuth matched filter, we can determine the target's moving direction. Take WAN HAI 506 for example. First, we compare the SAR images produced by two azimuth matched filters defined in (31) and (32).  (32), in which we can see that the target signal energy is concentrated together. Figure 16 indicates the SAR image with (31). Due to mismatch during the azimuth compression, the target signal energy in Figure 16 has been split into two parts. Therefore, the correct azimuth matched filter is (32). According to the derivation in Sections 2 and 3.5, (32) corresponds to (9), which means that the satellite and target are on different sides with respect to the receiver before the target crosses the LOS of the receiver antenna. Tables 2 and 3 provide the azimuth angles of the LOS of the receiver antenna and the satellite, respectively. From the azimuth angle relationship, we know that the satellite is on the left side with respect to the receiver, like in Figure 3b. Based on the above analysis, we can determine that the cargo ship is on the right side before crossing the receiver and moves towards the left side, which follows the direction of the arrow marked in Figure 12a.
Remote Sens. 2020, 12, x 19 of 21 above analysis, we can determine that the cargo ship is on the right side before crossing the receiver and moves towards the left side, which follows the direction of the arrow marked in Figure 12a.

Conclusions
In this paper, GNSS-SAR technique was proposed for maritime moving target imaging. The main work contained three parts. First, an approximate bistatic range history was derived to describe the phase variation of the target signal, the maximum error of which was less than 1 m compared with the real bistatic range. Second, a moving target imaging formation algorithm was proposed, involving unknown range cell migration correction, motion parameter estimation, and azimuth compression derivation. Finally, two groups of real experimental data were used to validate the effectiveness of the proposed method.
The obtained SAR images demonstrated that the target-to-receiver vertical ranges of the moving targets along the range had good accordance with the ground truth and the measurement results of the ships' length along the cross-range were very close to the real values, which confirmed the effectiveness of the proposed algorithm. Furthermore, Monte Carlo simulations showed that the MSEs of the proposed motion parameter estimation method were better than that of LSM in low SNR scenes. The capability of the proposed algorithm for the target's moving direction determination was also confirmed. In the case of a low sea state, target translation motion was dominant, while rotation motion could be negligible. The proposed method has the capacity to be applied for maritime moving target imaging, velocity estimation, and ship length measurement in coastal areas.
This paper considered a simple case of bistatic acquisition geometry where the receiver and the satellite are stationary. In the next work, more general cases should be taken into account. For example, when the receiver is mounted on the aircraft. In this case, Doppler ambiguity will occur. Moreover, the use of multiple GNSS satellites can further improve the information about target features. Thus, exploring data fusion of multiple SAR images should be investigated in the future.

Conclusions
In this paper, GNSS-SAR technique was proposed for maritime moving target imaging. The main work contained three parts. First, an approximate bistatic range history was derived to describe the phase variation of the target signal, the maximum error of which was less than 1 m compared with the real bistatic range. Second, a moving target imaging formation algorithm was proposed, involving unknown range cell migration correction, motion parameter estimation, and azimuth compression derivation. Finally, two groups of real experimental data were used to validate the effectiveness of the proposed method.
The obtained SAR images demonstrated that the target-to-receiver vertical ranges of the moving targets along the range had good accordance with the ground truth and the measurement results of the ships' length along the cross-range were very close to the real values, which confirmed the effectiveness of the proposed algorithm. Furthermore, Monte Carlo simulations showed that the MSEs of the proposed motion parameter estimation method were better than that of LSM in low SNR scenes. The capability of the proposed algorithm for the target's moving direction determination was also confirmed. In the case of a low sea state, target translation motion was dominant, while rotation motion could be negligible. The proposed method has the capacity to be applied for maritime moving target imaging, velocity estimation, and ship length measurement in coastal areas.
This paper considered a simple case of bistatic acquisition geometry where the receiver and the satellite are stationary. In the next work, more general cases should be taken into account. For example, when the receiver is mounted on the aircraft. In this case, Doppler ambiguity will occur. Moreover, the use of multiple GNSS satellites can further improve the information about target features. Thus, exploring data fusion of multiple SAR images should be investigated in the future.