Maritime Ship Target Imaging With GNSS-Based Passive Multistatic Radar

In the field of maritime surveillance, the global navigation satellite system (GNSS)-based passive radar has proven its potential for moving target detection (MTD), localization, and velocity estimation. The next stage is to investigate the possibility of obtaining the radar image of the moving ship for target recognition. However, the limited signal power budget of GNSS prevents the conventional inverse synthetic aperture radar technique that is based on target rotational motion and short observation time for GNSS-based passive radar imaging moving target. In this article, a two-stage imaging processing method relying on the target translational motion over a long observation time is proposed. The first stage confirms the presence of the target by a long-time MTD processing technique. In the second stage, based on the analysis of the Doppler history of the target signal in the slow-time domain, short-time Fourier transform and modified random sample consensus are combined to robustly estimate target velocity with reduced computation complexity. To obtain the focused bistatic image, azimuth compression is conducted by using the estimated target velocity. Finally, an image fusion operation is implemented to combine the bistatic images achievable from multiple satellites so that a multistatic image with high quality can be created. The effectiveness of the proposed method is confirmed by the real experimental results of three cargo ships illuminated by several satellites.


I. INTRODUCTION
S INCE the early 1990s, the signals of the global navigation satellite system (GNSS) reflected off the Earth's surface has been exploited to retrieve the geophysical and geometrical parameters of the sensed scene. This novel technology is the distinguished GNSS-Reflectometry (GNSS-R). Plenty of remote sensing applications based on GNSS-R have been developed, such as ocean altimetry [1], wind speed inversion [2], soil moisture retrieval [3], and deformation monitoring [4].
In recent years, the use of Earth-reflected GNSS signals for land or sea target detection has been investigated. However, due to the strong clutter and ambiguous target geolocation Manuscript  in the delay-Doppler map, the conventional GNSS-R with a forward-scattering configuration is unsuitable for such an application [5], [6]. A backscattering configuration is a measurement geometry preferable to the forward scattering configuration [5]. From a radar perspective, a system receiving the backscattered GNSS signals can be considered a GNSS-based passive radar system. As a passive radar system, it has some inherent virtues of low-cost, license-free, security, and no electromagnetic pollution, which provides a good complement to active radar systems. With the development of wireless communication technology, many terrestrial radio sources are employed as illuminators of opportunity for passive radars, such as digital video broadcasting-terrestrial (DVB-T), frequency modulated broadcasting, and Wi-Fi. Compared with terrestrial passive radars, GNSS-based passive radar has the advantage of persistent signal coverage over the world, particularly, in the coastal areas, open seas, and pole regions, thanks to the large constellations in space (global positioning system (GPS), GLONASS, Galileo, and Beidou). Moreover, the GNSS signals typically operate in the microwave region (L-band), which can easily penetrate clouds and fog. Therefore, compared with the conventional optical or infrared space-based sensors, GNSS-based passive radar is capable of working all day and all weather. Due to the abovementioned advantages, GNSS-based passive radar has attracted substantial attention from the radar research community. A series of relevant studies have been made over the last few years. One of the well-known GNSSbased passive radar techniques is the passive synthetic aperture radar (SAR), which concentrates on static objects on land, aiming to yield local area monitoring [7], [8], [9]. Moving to dynamic objects, in recent years, the application of maritime surveillance is a research hotspot. Two proof-ofconcept studies for maritime moving target detection (MTD) were implemented in [10] and [11]. Nevertheless, the main weakness of GNSS-based passive radar is the low signal power budget, restricting the radar's operational range. To address such an issue, some long-time integration techniques have been proposed for detection purposes [12], [13], [14], [15]. Because the technology of code (or frequency) division multiple access is applied in GNSS, an individual receiver (Rx) can separate the GNSS signals coming from different satellites. Noticeably, GNSS-based passive radar belongs to a multistatic radar in nature, enabling target location and velocity estimation. Relevant works can be referred to [16], [17], [18], [19], and [20]. The next step is to investigate the possibility of obtaining the radar images of the detected targets, which would bring an additional benefit, allowing target classification. For example, This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ the obtained images enable us to estimate the length of the ships and identify their characteristic features when different types of ships pass through the surveyed area. The inverse SAR (ISAR) technique applied in the active radar systems can be modified for the passive radar systems to realize moving target imaging. Some contributions to the ISAR images can be found concerning the exploitation of third-party illuminators [21], [22], [23], [24], [25], [26]. The capabilities of imaging naval and aerial targets by using DVB-T/DVBsatellite (DVB-S) signals have been shown in [21] and [22]. The experimental results of ground target imaging by using GSM and Wi-Fi transmissions have been reported in [23] and [24], respectively. For coastal surveillance, geostationary telecommunication satellites [25] and hitchhiking on the cooperative coastal radar [26] have been put forward to produce passive ISAR images, respectively. However, in the field of GNSS-based passive radar, a few studies have been presented to tackle the imaging issue of the moving target. Some preliminary experimental results of large-size ship imaging have been reported in [27] and [28]. Then, in [29], a proper processing chain to achieve the images has been developed. It can focus the energy of target signals in the bistatic range and Doppler domain to obtain the ship image. A mathematical framework is defined to map the bistatic range and Doppler domain into the equivalent monostatic range and cross-range domain for the interpretation of the image products. During the imaging processing, the proposed technique in [29] requires some accurate kinematic parameters of the ship target, such as the Doppler rate and the angle of the direction of arrival, retrieved by the stages of target detection and location. However, these parameters are usually achieved through a search process that is characterized by a heavy computational cost. Furthermore, most proposed algorithms of target detection and location assume that the ship behaves from an electromagnetic point of view as a point-like target whose reflection position is located at the "center of gravity" of the ship. On the contrary, in the real scenario, the reflection position on the target changing with the illumination and scattering angles may be randomly distributed within the target size during the observation time [17], [19], which would degrade the estimation accuracy of these kinematic parameters.
To image the moving target, the conventional ISAR technique uses the fast Fourier transform (FFT) in the azimuth direction for azimuth compression, which is valid for targets with pure rotational motion over a short observation time. However, due to the low signal power budget of GNSS signal, the limitation of the observation time is in contrast with the long aperture time needed by GNSS-based passive radar. On the other hand, the restricted signal power budget enables the major use of GNSS-based passive radar mostly suitable for short-range surveillance scenes, such as coastal areas and river channels, where the ship target moves with a translational motion and negligible rotational motion (i.e., low sea state conditions) [25]. Therefore, in these short-range application scenes, it is of interest to obtain the passive ISAR image relying on the translational motion over the long observation time, particularly, when the trajectory of the moving ship crosses the line of sight (LOS) of the radar antenna. This article is dedicated to the use of the GNSS signals for ship target imaging by resorting to the SAR technique because the ship translational motion is equivalent to the stationary Rx moving along the opposite direction with the same speed, similar to in strip-map SAR. Conventional GNSS-SAR algorithms [30], [31], [32], [33] have the knowledge of the radar platform velocity and then implement azimuth compression to fulfill imaging processing. However, in our case, the moving ship is a noncooperative target whose velocity is unknown. Hence, before azimuth compression, the target velocity must be estimated through the Doppler history of the target signal. In [34], motivated by the back-projection algorithm (BPA) commonly used in the field of GNSS-SAR, a moving target imaging method is proposed, in which the target signal energy focused in the bistatic range and Doppler domain is projected into the local Cartesian plane by searching the target velocity. However, like BPA, the projection process in [34] is performed at the expense of a high computational cost. In this article, to efficiently obtain the focused image, a two-stage imaging method is proposed. In the first stage, a long-time MTD processing technique is performed to determine the existence of the target in a range-Doppler (RD) map. In the second stage, a bistatic acquisition geometry of the imaging scene is first defined to analyze the Doppler history of the target signal with respect to the slow time. Then, based on the Doppler history of the target signal, short-time Fourier transform (STFT) [35] and random sample consensus (RANSAC) [36] are combined to estimate target velocity with reduced computation complexity and robust estimation errors. Finally, the azimuth-matched filter is designed by using the estimated target velocity to accomplish azimuth compression in the frequency domain. In our previous work [37], the preliminary images of ship targets were achieved via real experimental data. However, the inherent multistatic nature of GNSS-based passive radar is not taken into account, which can enrich image information from various illumination angles. Therefore, in this article, an image fusion operation is conducted to combine the bistatic images from multiple baselines, obtaining a multistatic image with enhanced quality. The proposed imaging method is tested against the experimental data, including several GPS satellites and three cargo ships with different lengths and speeds. The collected results validate the accuracy and robustness of the proposed target velocity estimation method and show that the multistatic image can improve the estimation accuracy of ship length.
This article is organized as follows. Section II describes GNSS-based passive radar detecting the moving target in the short-range surveillance scene. On this basis, Section III presents the processing procedure of velocity estimation and imaging formation. The field trial results are reported in Section IV. Finally, a conclusion is drawn in Section V.

II. LONG-TIME MTD PROCESSING
In this section, the 1-D GNSS signal model is first described. Then, the 2-D radar data formats of both the direct signal and the target signal are given, convenient for the radar signal processing. Finally, the long-time MTD processing technique is presented for the short-range surveillance scene.

A. Signal Model
Typically, the GNSS signal can be regarded as the product of three time waveforms [38]: the pseudorandom noise (PRN) code, the navigation code, and the radio frequency (RF) carrier. Taking the GPS L1 signal, for example, the simplified emitted signal can be expressed as where C(·) is the PRN code sequence whose duration is 1 ms with a chip rate of 1.023 MHz, D(·) is the navigation code, f c is the central carrier frequency, and ϕ ini is the initial carrier phase. Due to the noncooperative nature of the waveform of opportunity, the direct signal is usually employed as the local reference signal for the target signal compression in the passive radar system. Therefore, as shown in Fig. 1, the Rx has a pair of RF channels, comprising the reference channel (RC) and the surveillance channel (SC). The RC collects the direct signals emitted by the satellite (Tx) through a low-gain antenna, while the SC records the reflected signals from the target (Tg) moving on the surveyed area of interest through a high-gain antenna. However, the signal-to-noise ratio (SNR) of the direct signal output from the RC can be as low as −30 dB [39]. To overcome such a problem, signal synchronization is a crucial step that acquires the direct signal and tracks the observation parameters, including code delay, carrier phase, carrier frequency, navigation message, and so on. These parameters can generate a noise-free replica of the direct signal as the local reference signal. For the radar applications, the local reference signal stored as the 1-D data stream must be transformed into the 2-D radar data matrix in terms of the pulse repetition interval (PRI) that can be matched to the PRN code length. After intermediate frequency demodulation and ignoring constant phase and amplitude terms, the local reference signal can be written as where τ ∈ [0, PRI] is the fast time, u∈[−T /2, T /2] is the slow time, and T is the entire observation time.
τ di (u) = (R b (u))/(c), ϕ di (u) = 2π(R b (u))/(λ ), and f di (u) = −(Ṙ b (u))/(λ ) are the instantaneous time delay, phase, and Doppler frequency of the direct signal, respectively. R b is the baseline length between the Tx and the Rx, c is the speed of light, and λ is the carrier wavelength. τ e1 and ϕ e1 are the total delay and phase errors, respectively. They are induced by the atmospheric factors (i.e., troposphere and ionosphere delay) and the Rx artifacts (i.e., clock cycle slips and local oscillator drift). Likewise, supposing one point-like target, the target echo can be expressed as where τ r (u) = (R(u))/(c), ϕ r (u) = 2π(R(u))/(λ ), and f r (u) = −(Ṙ(u))/(λ ) are the instantaneous time delay, phase, and Doppler frequency of the target echo, respectively. R = R t + R r is the propagation range of the target signal, R t is the range between the Tx and Tg, and R r is the range between the Rx and Tg. τ e2 and ϕ e2 are the total delay and phase errors, respectively. It should be noted that the navigation code causes the phase transition, unfavorable to the target signal processing. Fortunately, the navigation codes in both (2) and (3) are the same within the range of 6000 km [15]. On the other hand, the delay and phase errors in both (2) and (3) are very close to each other due to the similar atmospheric factors and the shared oscillator between the RC and SC. The abovementioned interference factors can be eliminated later by range compression.

B. Long-Time MTD Processing Technique
The main goal of long-time MTD processing is to indicate the moving target in the RD domain and then extract the target signal with respect to the slow time through the constant false alarm rate (CFAR) technique [40]. The long-time MTD processing contains three main steps: range compression, range cell migration (RCM) correction, and long-time hybrid integration. Each step is described in the following.
1) Range Compression: Range compression can not only measure the time delay difference between the reference signal and the target signal but also eliminate the navigation code, the delay, and phase errors in the target signal. It is implemented in the fast time by cross-correlating the target signal in (3) with the complex conjugate of the reference signal in (2). After range compression, the time-domain range-compressed signal can be written as where CF(·) denotes the envelope of the cross correlation function and R bi = R t + R r − R b is defined as the bistatic range.
It should be noted that the calculation of (4) is conducted in the frequency domain based on the convolution theorem to realize high computational efficiency.
2) RCM Correction: As shown in (4), the envelope of the cross correlation function is coupled with the slow time. Consequently, the envelopes will be dispersed on several range cells, i.e., the RCM, with the ship translational motion during the long observation time. To align the envelopes within a single range cell resolution, the keystone transform (KT) is employed because it can correct the RCM without prior knowledge of target velocity. For detailed implementation steps, one can refer to [15]. After the RCM correction, the range-compressed signal in (4) can be rewritten as where R c bi is the constant bistatic range that the target appears to be after the RCM correction.
3) Long-Time Hybrid Integration: The energy gain offered by range compression is limited, which may be insufficient to extract the target signal from the background disturbance. Hence, to detect the target signal and reduce the computational complexity, a two-step long-time hybrid integration approach is performed here. In the first coherent integration step, the long observation time is divided into several frames (i.e., the range-compressed data matrix in (5) is divided into several small matrices with respect to the slow time). Each frame has an identical and short coherent processing interval (CPI), such as 2-3 s [13], [14]. Then, all frames conduct the FFT in a slow time. The target response in the lth RD map can be obtained as where r = τ × c is the bistatic range cell, l ∈ [0, N f −1], N f is the number of frames, f d ∈ [−PRF/2, PRF/2] is the Doppler frequency, PRF represents the pulse repetition frequency that equals 1/PRI, rect(u/T f ) is a rectangular window with a duration T f , T f is the length of CPI, W a ( f d − f l dc ) is the spectral envelop centered at f l dc , and θ l a (r, f d ) is the phase angle after the FFT. In a short-range surveillance scene, the coherent integration gain over a single CPI is sufficient to isolate the target response from the background disturbance because the moving ship yields a relatively high SNR input to the Rx [29]. In the second step, a noncoherent integration strategy is implemented to reduce the fluctuations of the background disturbance and increase the detection probability. The noncoherent integration of multiple frames is derived as As a result, a 2-D RD map is generated. It should be noted that the direct signal intruded from the sidelobes of the SC antenna in the real scene is the strongest return. After the long-time hybrid integration, the compressed direct signal (along with its sidelobes) is located at the zero bistatic range and zero Doppler frequency position. They can be filtered out by a simple direct signal cancellation procedure. Finally, the obtained RD map is used for CFAR detection. Fig. 2 shows the flowchart of the long-time MTD processing.  A simulation trial is performed to present the generated RD map, where a moving ship is modeled as five scatter points distributed uniformly along the target length and its trajectory is orthogonal to the LOS of the radar antenna. Additive white Gaussian noise is utilized as the background disturbance, while the intruded direct signal and the sea clutter are not considered here. The detailed simulation parameters are listed in Table I, where the GPS signal parameters are referred to [41]. Fig. 3(a) shows the resulting RD map. We can see that the track of the moving ship can be separated from the background noise and is aligned together at almost one bistatic range (R c bi = 2418m) since the RCM correction has been performed after range compression. For comparison, another example without the RCM correction is shown in Fig. 3(b), where the severe RCM can be seen in the bistatic range domain. Finally, in Fig. 3(a), the bistatic range of the track can be determined by the CFAR detection, which will be used in Section III-A.

III. MOVING SHIP IMAGING PROCESSING
In short-range surveillance scenes, such as coastal areas and river channels, the ship's translational motion is dominant in the case of a low sea state. The image of the moving target can be achieved by resorting to the SAR technique when the trajectory of the ship crosses the LOS of the radar antenna. In this section, an imaging method is proposed, mainly including the bistatic acquisition geometry analysis, target velocity estimation, and imaging formation. Each part is described as follows.

A. Bistatic Acquisition Geometry Analysis
As shown in Fig. 4(a), an east-north-up (ENU) reference system is constructed. The Rx is fixed at the origin, and the LOS of its radar antenna coincides with the due south direction. The Tx is in a state of quasi-static observing from the ground during the observation time (e.g., <120 s) because the distance between the Tx and the ground in an order of tens of thousands of kilometers leads to the Doppler bandwidth contributions caused by the target translational motion with respect to the Tx much smaller than the Doppler resolution cell (1/T f ) [29]. The trajectory of the Tg nearly orthogonal to the LOS of the radar antenna in the E-N plane is regarded as the synthetic aperture, which is equivalent to the Rx moving along the opposite direction with the same speed. Furthermore, the Tg is assumed to move with an almost constant velocity. This assumption is reasonable for the ship sailing at cruising speed [12], [13], [14], [15].
Considering the fact that R t and R b (>20 000 km) are much greater than R r (a few km), line Tx-Rx is almost parallel to the line Tx-Tg in Fig. 4. Therefore, the bistatic range defined in (4) can be approximated as where θ bi is the bistatic angle. From the geometry in Fig. 4(a), we know that θ bi has the trigonometric relationship as follows: cos θ bi = cos α Tx × cos ϕ Na−Tg−Rx (9) where α Tx is the satellite elevation angle and ϕ Na−Tg−Rx is the angle that varies with the target's motion. Through the top view in Fig. 4(b), the variation of (9) over the observation time can be further expressed as where u ∈ [0, 2L a /v], L a is the half-length of the synthetic aperture, v is the velocity of the Tg, ϕ Tx is the satellite azimuth angle (measured clockwise from the north), R s is the shortest range between the Rx and the Tg, and u a = L a /v is the slow time at which the Tg crosses the aperture center. It should be noted that, if the LOS of the radar antenna does not coincide with the due south direction, ϕ Tx should be recalculated. As shown in Fig. 5, it is not hard to understand that the recalculated ϕ Tx equals the absolute value of the difference between the azimuth angle of the LOS with the subtraction of 180 • and the real satellite azimuth angle (measured clockwise from the north).
Substituting (10) into (8), we get the bistatic range history over the observation time as It should be noted that, if the Tg in Fig. 4(b) moves along the opposite direction, the bistatic range history is changed as This peculiarity has been discussed in [37]. One can refer to [37] for detail. For the sake of convenience, only (11) is used for the next analysis. Expanding (11) into a second-order Taylor series approximation around u = u a , we have The Doppler history of the target signal can be derived by the first derivative of (13) with respect to u as and the Doppler rate γ can be derived by the second derivative of (13) with respect to u as Equation (15) indicates that the target signals from different satellites have a close Doppler rate value due to the satellites being in a state of quasi-static during the observation time. This peculiarity will be validated by the real experiments in Section IV.

B. Target Velocity Estimation
As evident from (14), the Doppler history of the target signal with respect to the slow time changes approximately linearly, i.e., a linear frequency modulation (LFM) signal. The target velocity can be obtained by estimating the Doppler rate value of the LFM signal. The long-time integration techniques reported in [12], [13], [14], and [15] are able to estimate the Doppler rate value. However, they are characterized by a heavy computational cost because of the search process. In this section, taking into account the computation complexity and accuracy, the STFT and the RANSAC are combined to estimate the target velocity.
Since the Tg has been detected at the bistatic range R c bi in the RD map, the target signal with respect to the slow time can be extracted from the range-compressed data matrix in (5). Then, the STFT projects the target signal into the time-frequency (T-F) domain by sliding a short-time window over the entire observation time. The STFT of the target signal with respect to the slow time can be written as where f w ∈ [−PRF/2, PRF/2] is the Doppler frequency, u ′ ∈ [−h/2, h/2], and win h (u) is the short-time window function with the length of h (i.e., the length of CPI). As shown in Fig. 6(a), a T-F representation simulated by using the parameters listed in Table I is achieved through the STFT with a window length of 2048 ms. We can see that the magnitude of the spectrum of the target signal distributed on the T-F plane forms an approximately straight line, which follows the LFM signal assumption. It is apparent that the slope of the straight line is the Doppler rate.
In the T-F representation, the conventional ways to estimate the Doppler rate are the Radon transform [38] and the Hough transform [42]. However, the computation burden of these methods is extremely heavy due to the LFM signal with a large time-bandwidth product. Linear fitting can be employed, which is simple and fast. The spectrum of the LFM signal should be first extracted by using the position of the STFT maxima [43] f w (u) = arg max Consequently, a T-F set S = {(u i , f w,i ) | i ∈ [1, N p ]} is obtained, as shown in Fig. 6(b). N p denotes the total number of points in the set. Because the target signal is contaminated by the background noise, the T-F set in Fig. 6(b) comprises both inliers, i.e., the points that can be fit to a line, and outliers, i.e., the noise points that cannot be fit to this line. A straightforward way to fit the line is the least squares method (LSM). However, the LSM optimally fits both the inliers and outliers, resulting in a bad fit to the T-F set. The RANSAC algorithm widely used in the field of computer vision is able to exclude the outliers and fit the line by using the inliers. Therefore, we modify the RANSAC algorithm suitable for the target velocity estimation. The proposed algorithm is referred to as the STFT-MRANSAC. It is summarized as follows.
Step 1: Implement the STFT via (16) and obtain the T-F set S via (17).
Step 2: Let K represent the number of iterations.
Step 3: Randomly select two points (u k 1 , f k w,1 ) and (u k 2 , f k w,2 ) from S to generate the line via the fitting model in (14). Superscript k denotes the kth iteration. Meanwhile, the Doppler rate γ k can as the candidate solution, and the Doppler centroid f k cen can be calculated by the above two points. To improve the robustness of the algorithm, γ k can is in a constraint given by where v min and v max are the minimum and maximum ship velocities, respectively. If γ k can does not satisfy the above constraint, the random selection process will be repeated.
Step 4: Calculate the distance of every point (u k i , f k w,i ) in S to the line generated in Step 3 by If d k i is less than the predefined tolerance threshold th tol , (u k i , f k w,i ) is considered an inlier. According to (14), th tol is defined as the Doppler difference between the bow and stern: where L s represents the ship length.
Step 5: Determine how many points from S fit the model within the given th tol . If the fraction of the number of inliers over the total number of points in S exceeds a predefined threshold th scale , γ k can is put into the candidate solution set C = {γ k can | k ∈ [1, K ]}. It should be noted that the abovementioned configuration parameters will be given in Section IV.
Step 6: Let k = k + 1. If k ≤ K , go back to Step 3. Otherwise, proceed to Step 7. Step 7: Use a median filter to determine the estimated Doppler rate value from C. Finally, the target velocity can be solved by (15). It should be noted that the estimation of R s will be explained in Section III-C.
The reason why a median filter is used can be interpreted as follows. The moving ship is no longer regarded as an individual point-like target in the short-range surveillance scene. Instead, it can be modeled as several scattering points distributed within the ship size. In other words, the target signal scattered from the moving ship is composed of multicomponent signals. During the translational motion, the spectra of multicomponent signals can form a wide line in the T-F representation owing to the Doppler frequency difference between the bow and stern. Fig. 7(a) shows the zoomed T-F representation whose Doppler resolution cell is about 0.49 Hz [highlighted by a red box in Fig. 6(a)]. We can see that the Doppler frequency difference between both edges of the wide line exceeds 5 Hz (i.e., more than ten Doppler frequency cells are shifted). However, the wide line brings difficulty to the fitting process. As shown in Fig. 7(b), it is possible to obtain an arbitrary Doppler rate that is located at the range between "Line A" and "Line B" by using the conventional RANSAC algorithm. In this case, the median filter solution (red line) can effectively increase the accuracy and robustness of the parameter estimation.
The STFT enjoys a low computational cost because it is performed on a computer by using the FFT with the calculation complexity of O(N log N ) operations. In Step 3, the random selection operation may be repeated one or M times. While the calculation of the distances between the points and the line in Step 4 is in the order of O(L). The operations in Steps 3 and 4 are conducted K times. Hence, the proposed STFT-MRANSAC has a maximum computational complexity of O[N log N + (L + M)K ]. Compared with the long integration techniques in [12], [13], [14], and [15] and Radon transform [38], where the complexities are in the order of O (N 2 log N ), the STFT-MRANSAC is more favorable from the calculation complexity point of view.

C. Imaging Formation
As the target velocity has been estimated in Section III-B, an azimuth compression operation can be performed to obtain the bistatic image of the ship. The bistatic image should be converted into an equivalent monostatic range and cross-range domain, convenient for the ship length estimation. Considering the inherent multistatic nature of GNSS-based passive radar, the bistatic images achievable from multiple satellites are combined to create a multistatic image. The abovementioned steps are described in the following.
1) Azimuth Compression: In the focused SAR processing, azimuth compression is fulfilled by removing the azimuth modulation of the target signal in the Doppler (or azimuth-frequency) domain with the matched filter. Therefore, designing the azimuth-matched filter is the key point here. An analytical expression of the azimuth spectrum of the target signal in the range and azimuth-frequency domain should be derived. For simplicity, we first consider an individual pointlike target. The azimuth FFT of the range-compressed signal in (5) can be written as (21) where f u ∈ [−PRF/2, PRF/2] is the azimuth-frequency. Substituting (11) into (21), we exploit the principle of stationary phase (PSP) to obtain an analytical solution on the integral in (21). After some lengthy algebra (see the Appendix), the azimuth spectrum is derived as (22), shown at the bottom of the page, where ρ af denotes the magnitude of the azimuth-frequency spectrum that is not important and will be omitted in the following analysis. As shown in Section III-B, all scattering points on the target surface experience the same Doppler history. Therefore, the azimuth spectra of multiple scattering points can be extended as (23), shown at the bottom of the page, where u n a represents the slow time at which the nth scattering point passes through the aperture center. From (23), it is understood that the modulation phase (the first term) should be removed. Hence, the azimuth-matched filter is the complex conjugate of the modulation phase as follows: Afterward, the azimuth compression process is done for each range cell by multiplying (23) with (24), followed by an azimuth inverse FFT (IFFT). The azimuth compression operation is given by As evident, the above implementation of azimuth compression is based on the FFT and IFFT. The total computational cost is in the order of O(N 2 log N ) for all range cells.
2) Range and Cross-Range Scaling: Range and cross-range scaling are required to transform the fast-and slow-time domains in (25) into the equivalent monostatic range and cross-range domain. From (13), we know that the bistatic range can be converted into an equivalent monostatic range if u = u a is chosen as the reference instant after the RCM correction. Since the Tg has been detected at the bistatic range R c bi , R s can be calculated by Correspondingly, the equivalent monostatic range is scaled by r x = (c×τ )/(1+cos α Tx cos ϕ Tx ). The cross-range is scaled by r y = v × u. Finally, the resulting image is obtained with an equivalent range and cross-range domain as follows: where R n a = v × u n a represents the cross-range position of the nth scattering points. Let us analyze the cross-range resolution.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. It is obvious that the maximum Doppler bandwidth (i.e., the Doppler history during the observation time) of the target signal is dependent on the antenna beam angle ϑ bw , described as Correspondingly, the optimal cross-range resolution is derived as and λ = 0.19 m. We can observe that the larger the antenna beam angle, the better the crossrange resolution. The cross-range resolution of a small antenna beam angle (<10 • ) is finer than the best range resolution (15 m) of Galileo E5a/b or GPS L5 signals, Therefore, the high cross-range resolution is sufficient to estimate the length of the maritime ship.
3) Image Fusion: As GNSS-based passive radar belongs to multistatic radar in nature, the bistatic images achievable from multiple baselines can be combined through image fusion to form a multistatic image. The multistatic image can acquire more scattering information about the moving ship due to the satellites illuminating the ship from various angles, favorable to the ship length estimation. The pixel-level fusion process is used here. However, image registration should be first performed to make the ship target in different bistatic acquisition geometries locate on the same pixel area. As shown in (25), the locations of the scattering points in the slow-time domain are only related to times when they passed through the LOS of the radar antenna. This implies that the homologous pixels in different bistatic images are almost located at the same pixel area in the cross-range domain, while the detected bistatic range bin cell values are affected by different radar configurations. Hence, image registration is performed to align the range cells concerning multiple baselines. Then, the noncoherent summation of multiple bistatic images is conducted, obtained as where m ∈ [1, M] and M is the number of satellites. Fig. 9 shows the overall moving target imaging and fusion process.

IV. EXPERIMENT
A maritime measurement campaign was undertaken at the Cyberport Waterfront Park in Hong Kong on May 16, 2019. Fig. 10(a) demonstrates the experimental site where the receiving hardware (red push-pin) was fixed on the shore, recording the direct signals from the satellites and the reflected signals from the moving ships, simultaneously. Fig. 10(b) shows the photograph of the receiving hardware that consists of the Rx front end, two different types of antennas, and the host computer. The Rx front end [see the inset in Fig. 10(b)] has two independent acquisition channels connecting to the circular antenna (CA) and the square antenna (SA). The CA is a commercial off-the-shelf right-hand circular polarization antenna for direct signal collection, while the SA is a custom left-hand circular polarization antenna with a high-gain in charge of the reflected signal collection. To validate the proposed method, three cargo ships with different lengths and speeds were selected as the targets of interest, whose trajectories were nearly perpendicular to the LOS of the SA. Fig. 10(c) illustrates the photographs of three cargo ships, including their names and the times when they passed through the field of view of the SA. The voyage-related parameters of three cargo ships provided by the automatic identification system (AIS) are listed in Table II, which can be used as the ground truth. The L1 band signals emitted by GPS satellites were exploited to obtain the bistatic images of three cargo ships. Table III gives the experimental and signal processing parameters. It should be noted that the former two cargo ships in Fig. 10(c), COSCO ROTTERDAM and KUO LIN, are chosen to give the detailed analyses, while, for WAN HAI 313, only the final image fusion result and the binary image will be shown to avoid description redundancy. Fig. 11 demonstrates the resulting RD maps of both COSCO ROTTERDAM and KUO LIN obtained from multiple baselines. It should be noted that the compressed direct signals (along with their sidelobes) as the strongest returns have been filtered out to enable the target responses visible in the RD maps. Due to the short range and large size, we can see from Fig. 11 that the long-time hybrid integration processing  TABLE II  VOYAGE-RELATED PARAMETERS OF THREE CARGO SHIPS   TABLE III  EXPERIMENTAL AND SIGNAL PROCESSING PARAMETERS can provide adequate integration gain for the target responses isolated from the background disturbance during the observation time although the Doppler frequency migrations are not corrected for further integration gain improvement. On the other hand, the Doppler frequency migrations present the tracks of the two cargo ships. They have been aligned together in the bistatic range dimensions. Furthermore, we can find that the signal energy backscattered from two cargo ships is much stronger for some specific scattering angles than the other perspectives. The main reason is that the fluctuations in the radar cross section (RCS) of the target are related to the illumination and scattering angles in the real scene.

A. Detection Results
Due to the tracks of two cargo ships aligned in the RD maps, their bistatic ranges with respect to multiple baselines can be measured at the maximum intensity values of the tracks after the CFAR detection. Consequently, the shortest range R s can be calculated by (26). Table IV lists the numerical results of the measured bistatic ranges and the estimated shortest ranges. We can see that the estimation errors between the estimated shortest ranges and the AIS ground truth given in Table II are much smaller than the optimal range resolution (∼150 m) of the GPS L1 signal. Furthermore, the differences among the estimated shortest ranges with respect to two cargo ships are less than the sampling resolution (18.3 m) of this signal processing system. Hence, after imaging processing, the focused targets in all bistatic images will be aligned in the monostatic range domains, convenient for image fusion.

B. Velocity Estimation Results
The operation of target velocity estimation can be performed after two cargo ships are detected in the RD maps. Fig. 12 shows the T-F representations of two cargo ships obtained from multiple baselines. Although the fluctuations in the target's RCS cause the magnitude of the spectrum of the target signal nonuniformly distributed on the T-F plane, we can  still see the approximate straight line with a negative slope value in each figure, following the LFM signal model. Taking COSCO ROTTERDAM as an example, we can observe from Fig. 12(a) and (b) that the Doppler histories in both figures are different due to their bistatic acquisition geometries. However, the slope values of both lines (i.e., the Doppler rate values) pertaining to sat. 1 and sat. 2 are very close to each other because the Doppler frequency contributed by the relative motion between the target and the Rx is dominant, while the Doppler frequencies induced by the target translation motion with respect to different satellites are extremely small and can be neglected. Similar results for KUO LIN can be seen in Fig. 12(c)-(e). The above findings from the real experimental data are well in line with the theoretical expectations [see (14) and (15)]. Moreover, apart from the spectra of the target signals, the T-F representations in Fig. 12 also contain the interference signals, such as the sea clutter (marked by the red boxes) and the sidelobes of the intruded direct signal (marked by the yellow boxes). The target velocity estimation accuracy will be affected by the interference signals when the conventional LSM is used for linear fitting.
The proposed STFT-MRANSAC is conducted to robustly estimate target velocity in the presence of many data outliers. First, the T-F set S is composed of the positions of the STFT maxima. Then, the MRANSAC can be performed. The configuration parameters used during the iteration of the MRANSAC are given as follows. The number of iterations K is set as 1000 so that a considerable number of candidate solutions (i.e., slope values) can be acquired. The minimum and maximum target velocities (v min and v max ) are 5 and 15 m/s, respectively, in order to have the candidate solutions in a constraint. The reason for the minimum target velocity configuration is that a too small value of the target velocity (such as 1 m/s) implies an increased observation time, during which the satellite may not satisfy the assumption of the state of quasi-static, while the configured maximum target velocity is high enough for most cargo ships moving on the sea surface. The tolerance threshold th tol for inlier determination can be calculated by the Doppler difference between the bow and stern. According to (20), L s is set as 300 m long enough for most cargo ships. Finally, the predefined threshold th scale for a candidate solution determination is set 5% to ensure that the candidate solution set C can be filled by as many candidate solutions as possible.
After the iterations for the fitting of straight lines, in order to better visualize the distribution of all candidate solutions, the clouds of points composed of the candidate solutions with respect to multiple satellites are illustrated in Fig. 13. It should be noted that, according to (15), the slope values representing the AIS ground truth are calculated via the parameters listed in Table II. We can see that the candidate solutions are randomly distributed within the constraint given by (18) due to the wide straight lines and the outliers. However, the AIS ground-truth data (black lines) are in the vicinity of the median filter solutions (red lines) in all the figures. This implies that the MRANSAC is capable of increasing the robustness of slope value estimation against various effects in the underlying interference signals and noise.    14 gives the fitting results of two cargo ships with respect to multiple baselines achieved from the median filter solution of the STFT+MRANSAC. For comparison, the generic RANSAC (whose code is from [44]) and the LSM are combined with the STFT to implement linear fitting. It should be noted that the configuration parameters of both the generic RANSAC and the MRANSAC are the same, and the constraint of the candidate solutions in (18) is also used in the generic RANSAC. We can observe that there are many outliers induced by the interference signals and noise in all the figures. As a result, the straight lines fit by the LSM (dashed blue lines) deviate significantly from the inliers due to the LSM optimally fitting all points. The straight lines fit by the generic RANSAC (green lines) are nearly located inside the inliers. However, the voting scheme adopted in [44] to find the optimal solution is not ideal in the wide straight line fitting issue. The estimation errors of the generic RANSAC may be worse than that of the LSM. On the contrary, the straight lines fit by the median filter solution of the STFT+MRANSAC (red lines) almost coincide with the lines generated by the AIS fitting lines (black lines) in all the figures.  Because the RANSAC estimates the parameters of a model by random sampling of observed data, ten trials are conducted to test the target velocity estimation performance. Fig. 15 depicts the target velocity estimation curves of two cargo ships with respect to multiple baselines provided by the abovementioned three methods. We can intuitively see that, compared with the STFT+RANSAC, the STFT+MRANSAC cannot only have reduced estimation errors with regard to the AIS velocities but also show a more robust estimation performance against the random selection procedure. The numerical results of the two methods are compared in Table V. For the STFT+LSM, its estimation performance is significantly affected by the outliers. Apparently, its estimation errors are larger than that of the proposed method in Fig. 15(b)-(e). However, there is one exception shown in Fig. 15(a). This can be explained from Fig. 14(a) where the dashed blue line is nearly parallel to the black line, but it only fits very few inliers. To obtain the final estimated target velocity for the

C. Imaging Results
As the target velocities are obtained, the passive ISAR images of the moving targets can be generated by using the proposed imaging method. Fig. 16 presents the resulting bistatic images of two cargo ships. The color scales are in dB, where 0 dB is the mean power of the noise. In all the figures, the responses of two targets can be distinguished from the background, and their positions are aligned in the range and cross-range domains. Comparing Fig. 16(a) with (b), we can see that both of them observe different components of the ship body due to the different bistatic angles. Similar results can be seen in Fig. 16(d) and (e). If the target features are extracted from these bistatic images for the ship length estimation, it is obvious that the estimated lengths are far from the actual ones. On the other hand, in Fig. 16(c), the low level of SNR provided by sat. 3 makes the outline of the ship target blurred. It seems that sat. 3 is less suitable for the imaging task. The above findings in the bistatic images indicate that it is essential to increase the information gathered about the ship target by combining the bistatic images achievable from multiple illumination angles.
Image fusion operation is conducted to combine the bistatic images for multistatic image generation. Fig. 17 shows the meaningful multistatic images of COSCO ROTTERDAM, KUO LIN, and WAN HAI 313. Compared with the bistatic images, the multistatic images have lower fluctuations of the background, making the outlines of the targets more visible. The ship length can be measured through the distance between both edges of the ship body along the cross-range direction. We can intuitively see that the ship lengths in Fig. 17(a)-(c) are, respectively, more than 250 m, less than 200 m, and a little more than 200 m, which are will in line with the ship length comparison listed in Table II. This reveals that the multistatic images allow roughly estimating the ship length and identifying the ship features for classification.
The binary images of three cargo ships are created to evaluate the numerical results of the ship lengths. According to the color scales in Fig. 17, a simple threshold is set as 10 dB to extract the target features. Fig. 18 shows the binary results. The ship lengths can be measured as the distance between the rightmost and leftmost pixels of the main ship body. While some isolated pixels out of the main ship body are discarded. The measured lengths of COSCO ROTTERDAM, KUO   two bistatic images are used for image fusion, resulting in low-intensity pixels belonging to the target, difficult to be isolated from the surrounding area. Consequently, the pixels not belonging to the target have more chances to be extracted by the threshold and are used for the length estimation. This problem can be overcome by exploiting more opportunistic signal sources. Moreover, the threshold selection has an effect on ship length measurement. Due to the weaker responses from the ships' bows and sterns shown in Fig. 17, the higher the threshold value, the shorter the measured length. Therefore, how to extract the target feature for accurate ship length estimation may be investigated in the future.
Taking KUO LIN for example, a simple sensitivity analysis is provided, which evaluates the imaging performance degradation against errors on the estimated target velocity. It should be noted that the imaging performance is quantitatively calculated via the image entropy in [34]. The estimated target velocity is defined asv = v 0 + e v , where v 0 is the true value of the target velocity in Table II, and e v is a zero mean Gaussian random variable with a standard deviation   Table VI shows the standard deviation (accuracy) of the image entropy error as obtained by 10 3 independent realizations ofv. We can see that, in all cases, the standard deviation values increase with the values of η v increasing. Therefore, the imaging performance relies on accurate target velocity estimation.
Because some imaging methods tailored for GNSS-based passive radar have been published in recent articles [29], [34], it is necessary to discuss the advantage and disadvantages of the proposed method. As shown in Sections II and III, the most complicated operations in the proposed method are the RCM correction, velocity estimation, and azimuth compression, whose total calculation complexity is in the order of O(N 2 log N ); for the existing method in [29], if the calculation cost of the motion parameters used for the scaling steps can be ignored (i.e., they are provided by the AIS), a cost in the order of O(N 3 log N ) is needed to focus the target energy in the RD map. While the cost of the projection process in [34] is more than O(N 3 ) but less than O(N 4 ) since the particle swarm optimization algorithm is adopted to reduce the computational burden. From the above analysis, we can see that the computational cost of the proposed method is lower than the existing methods in [29] and [34]. Due to the projection process, the image resolution in [34] is restricted by the RD map. Hence, taking KUO LIN for example, Fig. 19 shows the multistatic image via the existing method in [29] by using the experimental parameters given in Table III. Compared with the result in Fig. 17(b), it seems that the ship body in Fig. 19 shows a higher SNR. From Fig. 8 and the beam angle of the SA listed in Table III, we know that the cross-range resolution of the proposed method is relatively high (<1 m), while the cross-range resolution in [29] is dependent on the length of CPI. According to (17) and (18) given by Pastina et al. [29], the cross-range resolution in Fig. 19 is calculated as 9.3 m. Therefore, the high cross-range resolution of the proposed method causes the target energy to be dispersed in several cells. As a result, the ship body in Fig. 17(b) shows a lower SNR than that in Fig. 19. In summary, the proposed method has a low computational cost and high cross-range resolution. However, it is based on the requirements that the moving ship provides a relatively high SNR input to the Rx and the trajectory of the moving ship nearly perpendicular to the LOS of the Rx antenna. A river shipping monitoring scene may fulfill both requirements.

V. CONCLUSION
This article exploits the GNSS signals to obtain the passive ISAR image of the ship relying on the target translational motion over the long observation time in the short-range application scene. Therefore, a two-stage imaging processing method is proposed. In the first stage, a long-time MTD processing technique is conducted to determine the presence of the target in the RD map. In the second stage, the bistatic acquisition geometry of the imaging scene is analyzed to describe the Doppler history of the target signal induced by the translational motion during the long observation time.
Then, based on the Doppler history, the STFT+MRANSAC is presented to robustly estimate the target velocity. With the estimated target velocity, azimuth compression can be fulfilled by the designed azimuth-matched filter so that the focused bistatic image can be created. To obtain the multistatic image with improved quality, an image fusion operation is performed. Finally, the effectiveness of the proposed method is validated by the collected experimental data of two cargo ships illuminated by several satellites. The velocity estimation results indicate that, compared with the STFT+RANSAC and STFT+LSM, the proposed STFT+MRANSAC can enhance the accuracy and robustness of velocity estimation against various effects in the underlying interference signals and noise. The multistatic images show that it is possible to classify different types of cargo ships passing through the surveyed area in terms of the target length estimation.
The proposed imaging method requires the trajectory of the moving ship nearly perpendicular to the LOS of the Rx antenna. This requirement can be reached concerning river shipping monitoring but may not be easy on the sea surface. Therefore, the next work will consider more general bistatic acquisition geometry. Moreover, the ship length estimation can provide preliminary target identification. More target features need to be studied through a large number of experimental data. Therefore, more experiments will be carried out in the future.

APPENDIX
The PSP is employed to derive an approximation of the analytical representation of the integral in (21). To apply the PSP, a general case is first given as where A(u) is the slowly varying envelope function and ψ(u) is the phase term function. The main idea of the PSP is that the integral values of (A1) on both sides of the stationary point (the point where the derivative of the phase is zero) are canceled out due to the slowly varying or constant magnitude of the envelope. Hence, we should find the stationary point. This is obtained by computing the derivative of ψ(u) with respect to u and setting it to zero, i.e., by solving the following: where u * denotes the stationary point. Substituting the solution of u * to (A1), we achieve where C is a constant amplitude term. C and the constant phase offset π/4 can be omitted in the following derivation. Moving back to (21), the above procedure will be applied. Substituting (11) into (21), the stationary point u * is derived as Substituting (A5) back to (21) with rearrangement, we achieve the final analytical solution given by (22).