Maritime targets velocity estimation in space-based passive multistatic radar using long integration times

Target detection by means of space-based passive radar sensors generally requires the adoption of long integration time strategies to reinforce sufficiently the signal strength. These are usually based on the recovery of the target Doppler-rate to cope with the range and Doppler migration experienced over the long dwell considered. In this work, we put forward a Taylor-series approach that capitalizes on the set of bistatic Doppler-rates estimated in Global Navigation Satellite Systems (GNSS)-based multistatic radar configurations to estimate the velocity of ship targets with increased accuracy with respect to conventional Doppler-based procedures. Both the cases of single-element and array receiver configurations have been considered. Theoretical and numerical results under different use cases show as leveraging on the long integration times adopted at the detection stage could significantly increase the accuracy of the estimated ship velocity components. Few experimental results are also provided, verifying the potentialities of the proposed approach in operative scenarios of practical interest for this technology. The proposed approach is not limited at the GNSS case, but it could be potentially applied to any multistatic passive radar system.


I. INTRODUCTION
Passive radar sensors relying on satellite illuminators are a promising option for maritime surveillance purposes [1]- [4]. As passive systems, they do not rely on a dedicated transmitter, but they exploit RF signals already existing, thus allowing for sensors with reduced size, weight and complexity. Additionally, illumination from the space can assure coverage even in remote areas where terrestrial signals are absent (e.g., in open sea). Furthermore, constellations of satellites make possible operating with multistatic radar configurations by using a single receiving station [5]- [7].
The classical approach for passive radar detection consists in looking for signal peaks in the differential delay (range) and Doppler plane after cross-correlating the reference signal (often recorded by an additional dedicated channel) with the reflected signals from the surveyed area. However, one fundamental issue of space-based passive radar systems is that the power budget may be very restricted. The general trend to overcome this issue is extending the integration time, in order to strengthen the signal power sufficiently to enhance the target detectability [8]- [10]. Such class of methods consider additional search variables other than the conventional bistatic range and Doppler to cope with the target migration experienced over the integration time. Therefore, target detection occurs in a higher dimensional space, providing additional radar observables. Although these have been typically included in the detection space only to increase the power budget, they represent an additional source of information that could be exploited to enable advanced radar tasks. In this work, we investigate as additional radar observables provided by long integration time approaches can be fruitfully exploited in space-based radar systems for ship targets velocity estimation. The specific focus is on Global Navigation Satellite Systems (GNSS)-based multistatic radars, but the approach is suitable for different illuminators.
Over the last few years, several methods for extending considerably the integration time in GNSS-based passive radar have been proposed. These could be fully coherent [11] [12] or hybrid (i.e., multi-frame) [8] [13] [14]. In both the cases, most of the approaches implement a target motion compensation (TMC) procedure driven by hypothesized Doppler-rates: such information is therefore provided as detection output along with the bistatic range and Doppler shifts.
Detection outputs from multiple bistatic pairs are usually combined to determine the moving target location and velocity. Considerable efforts have been carried out to explore as multi-bistatic range and/or Doppler measurements can be fused to estimate the target location [5]- [7], [15], [16], while much less works explored the velocity estimation task. Methods for the joint estimation of the target position and velocity have been investigated in [18]- [19]. In [20], velocity estimation is accomplished by using multiple Doppler shifts in combination with the target position determined via a neural network algorithm in OFDM-based passive radar networks. In [21], a sparse Bayesian learning approach for the motion estimation of multiple targets is proposed. With specific regard to GNSS-based passive radar systems, velocity estimation via multistatic Doppler measurements is discussed in [22], while in [23] it has been experimentally verified that the velocity components of a ship target can be estimated by solving a linear system of equations relating the target velocity vector to the set of Doppler shifts.
All the above-mentioned velocity estimation methods rely on the different bistatic Doppler shifts estimated from the multiple transmitter-receiver pairs. Nevertheless, such Doppler-based approaches may result in degraded performance when the target is detected on a limited number of bistatic links, particularly when these links do not provide sufficient or suitable angular diversity. However, an accurate recovery of the ships velocity is a fundamental task for a reliable maritime situational awareness, which allows taking timely reactions against possible threats. Moreover, advanced system functionalities, such as imaging procedures, require a very accurate knowledge of the ship dynamics for recovering essential information for classification procedures [24]. Therefore, it is of interest to explore methods that capitalize on the Doppler-rates provided by the long-time approaches to improve the accuracy of the velocity estimation task.
To this purpose, a technique jointly using Doppler frequencies and Doppler-rates and based on Taylor expansion is here proposed. Some results along this line were anticipated in [25]. Extending such preliminary work, here we consider two possible configurations for the receiver: a first one operating with a single receiving antenna and a second one operating with an array of antennas. Such receiver configuration based on an array has been found adequate for the passive radar thanks to the possibility to provide angle of arrival measurements, usually unavailable due to the need to operate with staring antennas configurations to enable long integration times procedures [26] [27]. Furthermore, a weighted solution of the estimation problem is here considered, in lieu of the non-optimized combination of the measurements information previously developed. This is particularly relevant for space-based multistatic radar configurations, as the accuracy of the observables provided by different baselines can significantly vary due to several reasons, in particular the unpredictable variation of the bistatic radar cross section under widely separated illumination angles [6][7] [24].
In order to derive a theoretical benchmark for the performance, the Cramér-Rao Lower Bound (CRLB) is derived. Simulated results verify as the proposed technique provides performance close to this bound and may outperform the conventional Doppler shift based technique under different realistic situations. Moreover, experimental results comprising Galileo satellites as illuminators and targets belonging to different classes are reported, showing the effectiveness of the approach in applications of practical interest.
Therefore, the main contributions of this work can be summarized as follows:  A novel approach for ship velocity estimation in spacebased multistatic passive radar systems is proposed, which capitalizes on the Doppler-rate measurements made available by the long-integration time procedures adopted at the detection stage;  the theoretical performance of the method have been provided in closed-form via the Cramér-Rao Lower Bound, and compared to the performance achievable via conventional Doppler-based procedures;  the performance of the conventional and proposed approaches has been analyzed via simulations for a simpler receiver configuration using a single element radar antenna and a more efficient receiver using and array antenna under different use cases;  in order to prove the concept and verifying its practical potentialities, experimental results concerning river shipping monitoring and port operation applications have been provided. The rest of the paper is organized as follows: Section 2 provides a geometry description and recalls the conventional techniques relate to single frame processing, while Section 3 details the velocity estimation technique capitalizing on long integration times. Section 4 proposes a theoretical performance analysis through the CRLB, and results obtained against synthetic data are presented in Section 5. Experimental results are then shown in Section 6 and Section 7 closes the paper.

II. GNSS-BASED MULTISTATIC RADAR SYSTEM
The operative conditions comprise navigation satellites as illuminators of opportunity and a receive-only device in a remote location above the sea; for example, it could be mounted on the coast or, as depicted in Fig. 1, on a moored buoy. GNSS constellations guarantee that multiple satellites are simultaneously in visibility: as they operate on multiple access schemes, the receiver can easily discriminate the signals pertaining the different satellites, thus inherently forming a multistatic radar system. The receiver is equipped with two RF channels, one for recording the direct signals emitted by the satellites (reference channel) and another collecting the reflected signals from the surveyed area (radar channel). Usually, the radar channel employs a staring high-gain antenna, as rotating antenna beams are not adequate for the passive radar use due to the consequently limited time-on-target, contrasting with the requirement of sufficiently long coherent processing intervals (CPIs). System performance can be enhanced by adopting an array of antennas, which can allow increasing the available signal-to-noise ratio (SNR) as well as providing an estimate of the target direction-of-arrival (DoA). Both the cases of radar antenna using a single element and an array configuration will be considered in the following. and denote its elevation angle (measured with respect to the receiver position) and aspect angle (measured clockwise from -axis). Without loss of generality, the receiver position coincides with the origin and the -axis defines the radar antenna boresight.
In order to provide maritime situational awareness, the sensor has to detect the ship targets and, in turn, to estimate their position and velocity. The conventional approach to achieve these goals is depicted in Fig. 2. Range compression is implemented by cross-correlating radar and reference channels data (due to the low power budget, in the GNSS case a noise-free replica of the direct signal is usually reconstructed via signal synchronization [3]). Then Moving Target Indication (MTI) is performed via Doppler filtering followed by a detection stage extracting peaks in the resulting range and Doppler (RD) map. Thus, for each satellite an estimate of the bistatic range and Doppler shift of the target can be provided.
Particularly, for the system under consideration, the i-th bistatic range is equal to while the i-th Doppler shift is given by where is the wavelength (for sake of simplicity, here we assume the N signals sharing the same carrier frequency, however the approach can be immediately extended to the multi-frequency case). Then, an estimate of the target position can be obtained by intersecting the bistatic ranges [5]. In the case of an array receiver configuration, localization accuracy can be also improved by intersecting the N bistatic range ellipses with the target bearing information [26]. This is particularly useful in the case of poor localization accuracy provided by the bistatic ellipses intersection due to unfavorable geometric configurations. (In the case of multiple targets, a data association among the different bistatic links is in order [21], [28]. This task is beyond the scope of this work and hereinafter we will assume it successfully achieved).
Lastly, the velocity estimation task can be accomplished by exploiting the linear relation between the target velocity vector components and the bistatic Doppler shift. The Doppler-based approach presented in [23] is here briefly recalled. The bistatic Doppler (2) can be rearranged as  An estimate of the target velocity, ̂, can be therefore obtained by replacing in (4) , , and with their noisy values available after MTI and position estimation stages and then inverting the system (pseudo-inversion is implemented if N is such that the system is overdetermined).
In the remainder of the work, we will take the velocity estimation method in [23] described above as representative of the conventional approaches relaying on the Doppler shifts.

III. VELOCITY ESTIMATION TECHNIQUE USING LONG INTEGRATION TIME
It is worth noticing that the target kinematic parameters estimation stages in Fig. 2 implicitly assume the success of the detection task. However, the limited power budget of space-based passive radar systems makes the conventional MTI approach described above effective only for ships with sufficiently high radar cross section and/or at short ranges. Long integration time strategies are an effective solution to address the detection of less reflective ships and/or at higher ranges. As the processing interval increases, received signal may exhibit range walk and, for sufficiently wide bands, even range curvature. Moreover, the Doppler can be no longer assumed as constant but it might be characterized by a non-negligible Doppler-rate. In the geometry under consideration, for the i-th baseline this is given by Therefore, advanced MTI processors must be implemented to handle the range and Doppler migration, which is in the scope of target motion compensation (TMC) procedures. A long-time hybrid integration technique has been introduced in [13] that aims at accumulating the target energy over RD maps pertaining M frames (CPIs) by compensating the migration inside each frame and among consecutive frames. To this purpose, as shown in the top-part of Fig. 3, the whole integration time is segmented in consecutive windows, each with duration such that the target complex reflectivity can be assumed as constant. Then, TMC is implemented by using a bank of filters, each matched to an admissible value of the Doppler-rate, to correct the intra-frame and frame-toframe range and Doppler migration. The motion-compensated maps are then quadratically integrated, thus achieving a longtime RD map for each hypothesized Doppler-rate. The map corresponding to the ̇ value closest to the actual value provides the highest gain. Therefore, for each satellite the target can be detected in a data cube, thus providing a measure of its range, Doppler and Doppler-rate.
In this work, we aim at enhancing the velocity estimation task by capitalizing on such extra information. We refer to this method as Doppler+Doppler-rate velocity estimation technique, as both the sets of Doppler and Doppler-rate measurements are exploited to reach the objective. The block diagram of the method is depicted in the bottom part of Fig. 3.
Let ̂ be a rough estimate of in its vicinity; (5) can be approximated in Taylor series as follows where ∇ is the vector differential operator. The information space of the velocity estimation problem (4) can be augmented with the N pseudo-linear equations (6). After simple manipulations, we get Solving (7) can provide the shift Δ from the guess ̂ and the actual velocity vector. However, as (6) are pseudo-linear, an iterative procedure has to be implemented to reach the convergence. Let k be the index denoting the iteration number and let the apex (k) denote a variable value at the k-th iteration. Let ̂( −1) be the velocity guess available at the previous iteration, from which ( −1) and ( −1) are built. Particularly, the estimate of the position ̂ and ̂( −1) are used to set ( −1) , while Δ̇( −1) entry of ( −1) is achieved by correcting at each iteration the Doppler-rate measure obtained at the detection stage with the Doppler-rate evaluated in correspondence of ̂( −1) (i.e., (5) with =̂( −1) ). By pseudo-inversion, the current value of the estimated velocity shift Δ̂( ) is obtained and it is used to update the velocity estimate, i.e., ̂( ) =̂( −1) + Δ̂( ) . The system is initialized with the velocity estimate obtained with (4), i.e., ̂( 0) =̂, and it is iterated until the new correction becomes smaller than a certain threshold (or k equals a maximum value). We denote the velocity estimate at the output of this procedure as ̂,.
It is worth to point out that (7) is a mixed-mode system, since it combines measurements of Doppler (therefore, Hz) and Doppler-rate (Hz/s). The adoption of a Taylor-series estimation method is suitable for this type of problems, as this class of methods straightforwardly handles the combination of measurements of different types [29]. Particularly, in the case under consideration, the Doppler measurements contribute to the initialization of the algorithm, while the set of Dopplerrates regulates the refinement. Moreover, the Doppler information plays a fundamental role in the iterative procedure by means of the first N rows of the pseudo-linear system (4), i.e., Δ = . This constraints the subsequent refinements of the velocity estimate to lie on a particular direction according to the Doppler planes intersection due to the particular geometric configuration, thus improving the robustness of the approach against inaccurate measurements.
It is also worth to stress that the method combines independent measurements pertaining multiple bistatic acquisitions. Each measure is affected by its own error , usually modelled as a zero-mean Gaussian random variable, i.e., ~(0, 2 ), being 2 the variance of the distribution [7]. As it will be better detailed in Section 4, different system parameters affect the accuracy of the Doppler and Dopplerrate measurements. Moreover, the accuracy of a particular observable can vary with the particular bistatic link. This is particularly relevant in space-based passive radar systems, as different satellites can be characterized by widely spaced illumination angles, which could entail strong variations of the e.m. response of a ship target and ultimately of the power levels provided at the receiver. Therefore, the combination of the independent information should be optimized according to the expected accuracies of the exploited information.
To this purpose, Δ̂( ) can be obtained as the least-square solution of (7) by weighting the terms during the system inversion according to the covariances of the measurement errors [29]. Let Σ be the error covariance matrix of the Doppler and Doppler-rate measurements, at the k-th iteration the maximum likelihood velocity estimation correction is computed as So far, the velocity estimation approaches have been detailed for the general case of targets undergoing a threedimensional motion. However, as the focus of this work is on ship targets, we can assume the motion occurring in a 2D domain, i.e., component of is constant and therefore, = 0. Without loss of generality, we can assume = 0. In such conditions, the dimensionality of the estimation problems (4) and (7) can be reduced from 3 to 2 by deleting the third column of the system matrices and , respectively, and considering in (7) ∆ = [Δ ∆ ] . For clarity of illustration, we will focus on such simplified geometry model in the remainder of the paper.

A. CRLB DERIVATION
In order to derive a benchmark for the performance of the proposed method, the Cramér-Rao Lower Bound (CRLB) is here derived. Specifically, by referring to the scheme in ( | ) = 1  ̂ (resp. ̇) modelled as a N-dimensional Gaussian random variable with mean value (resp. ̇) and covariance matrix Σ = diag{ 1 2 , … , 2 } (resp. Σ̇= diag {̇1 2 , … ,̇ 2 }); the joint pdf of the available measurements can be written as in (9) being = 0 for the conventional Doppler-based method and = 1 for the proposed approach. The derivation is reported in Appendix. The CRLB is finally obtained by inverting . By means of the Woodbury formula, can be written as where , = , −1 and ̇= ( + ,̇) −1 ,̇, , being the conformable identiy matrix.
From (12) we observe as the bound is composed by two terms. The former is the CRLB concerning the conventional approach exploiting Doppler information only to recover the target velocity (provided that , is an invertible matrix), while the latter takes into account the impact on the performance of the exploitation of the Doppler-rate measurements (i.e., of the performance improvement as it will be shown in the next sub-section).

B. CRLB ANALYSIS
The comparison between the Doppler and Doppler+Dopplerrate based techniques is here carried out by referring to the (22) sub-matrix of related to the velocity components and evaluating the area of the confidence ellipse (i.e., the area of the iso-contour of the 2D Gaussian distribution) associated to a probability value equal to 0.99. The analysis is done by assuming:  the covariance matrix Σ (of the estimated position) equal to the CRLB evaluated using range measurements only (̂ with i=1,…,N) or by combining range and DoA measurements ( and with i=1,…,N). These measurements are supposed to follow independent Gaussian pdfs, namely ~( 0 , 2 ) and ~( 0 , 2 ), with standard deviations respectively equal to, [30], [27]: where (i) c and B represent, respectively, the speed of light and the satellites bandwidth; (ii) K, d and θ 0 are respectively the array elements, the inter-element distance, and the actual target DoA; and, finally, (iii) SNR i is the final signal-to-noise ratio reached at each bistatic link evaluated as SNR i = K√MSNR f with SNR f single-frame SNR (obtained over a CPI T), and M the number of RD maps non-coherently integrated to achieve the long RD map over which detection occurs.
 Doppler and Doppler-rate standard deviations given by, [31], [13]: where ̇ is directly derived from the CRLB for the Doppler rate, [31], considering an equivalent time interval equal to √M and is directly evaluated from ̇ as a consequence of the Doppler migration correction performed in the long time integration driven by the Doppler rate value, [13].
The CRLB analysis is here done by referring to the Galileo constellation and considering E5a bandwidth (carrier frequency: 1176.45 MHz, bandwidth: 10.23 MHz): firstly = 2 satellites are assumed available, namely the minimum number required to estimate the target kinematic state in case of ship targets. Both the single antenna and the array receiver configurations are considered. In the latter case, we set = 14 elements with inter-element distance equal to = 2 = 0.5 m: this ensures grating lobes at ±30° of steering direction, which is supposed to be sufficiently larger than the antenna beamwidth [27]. The ship target is assumed located at 0 = [3000,100] m, moving with velocity 0 = [5,10] m/s and detected with the long integration time approach by noncoherently integrating M = 10 RD maps ( = 3 s). The same SNRi (15 dB) is assumed at both the bistatic links and, to allow a fair comparison, for both the receiver configurations. Fig. 4 shows the ellipse area associated to the CRLB (in dB) as a function of the azimuth (∆ψ) and elevation (∆ϕ) diversity of the second satellite with respect to the first one, having fixed this last in [ϕ 1 , ψ 1 ] = [45°, 45°]. Particularly, Fig. 4 (a) and Fig. 4 (b) refer respectively to the Doppler-based and the Doppler+Doppler-rate based technique for the single antenna receiver, while Fig. 4 (c) and Fig. 4 (d) concern the array receiver configuration. By comparing the performance of the Doppler+Doppler-rate based technique to that of the Dopplerbased [ Fig. 4 (b) to Fig. 4 (a) and Fig. 4 (d) to Fig. 4 (c)] it is easy to observe that the exploitation of the Doppler-rate information provides a drastic reduction in terms of CRLB area corresponding to a better velocity estimate. Noticeably, the exploitation of the Doppler-rate measurements enables the estimation of the target velocity also in situations where the Doppler-based approach fails due to a not suitable angular diversity. This is easily observed by noticing that the "U" shaped area in Fig. 4 (a) (corresponding to the locus of the points where , is a singular matrix) is reduced to a narrow spot area around ∆ψ = ∆ϕ = 0 in Fig. 4 (b): such residual spot area is basically due to the fact that, as the angular diversity between the two satellites reduces to 0°, it is not possible to localize the target on the basis of range measurements only. Therefore, this area disappears when the array antenna receiver is used, Fig. 4 (d), as the availability of the DoA measurements allows to overcome this problem. For this last configuration, we observe that very good performance can be achieved on the whole (∆ , ∆ ) plane.
In order to analyze the influence of the value of Dopplerrate in the estimation accuracy, we considered the ship in the same position of the previous case study sailing at the same speed (‖ ‖ = 11.18 m/s) but varying its heading [defined as ℎ = 1 ( ⁄ )] between 0° and 90°. The different headings correspond to different values of the velocity components (i.e., different motion conditions and Doppler-rate values). The violet '□' markers in Fig. 5 show the ellipse area associated to the CRLB as a function of the Doppler-rates achieved for the different headings selecting = 2 satellites with positions  It is worth to recall that the main highlight of GNSS is the large number of available satellites, potentially offering the chance for selecting the more appropriate transmitters for addressing specific radar tasks, such as the velocity estimation considered in this work. Fig. 6 shows the ellipse area associated to the CRLB as a function of the number of selected satellites, being the first satellite located in [ 1 , 1 ] = [45°, 10°] and the others at = 1 and = 1 + ( − 1)∆ with ∆ = 10° and = 2, … , . For this analysis, a target located quite close to the receiver has been considered, 0 = [600,100] m, as this represents a very challenging condition. From the figure, it can be observed as the exploitation of the Doppler-rate brings improvements in the estimation accuracy for both the receiver configurations also when using more than two satellites, even if this improvement decreases as the number of satellites increases. This is due to the fact that the exploitation of more satellites provides a higher spatial diversity which gives rise to a quite accurate velocity estimation even by using the conventional (Dopplerbased) approach. However, this analysis shows that the improvement provided by the Doppler-rate in the velocity estimation compared to the conventional method is remarkable especially for cases of practical interest, namely when the target has been detected on a few bistatic links ( ≤ 4).

V. PERFOMANCE ANALYSIS
In this section, the effectiveness of the proposed approach is analysed by means of simulated data. All the shown results have been obtained by considering as a study case a ship target with size roughly equal to 40 m by 20 m sailing with velocity = [5,10] m/s. . This case can be assumed as a benchmark for the achievable performance of both the Doppler/Doppler+Doppler-rate based velocity estimation techniques.  Case B: as case A, but with estimated position. In this case, the performance of the localization will influence the estimation of the velocity since the knowledge of the target position is required by both the Doppler and Doppler+Doppler-rate based approaches.  Case C takes into account the possible variation of the target e.m. response with the illumination angle [6], [24]. To this purpose, the ship is assumed again to behave as a point-like target whose actual position, at each bistatic link, is modelled as a random variable uniformly distributed within the target size and uncorrelated among the baselines. As for case B, the velocity estimation performance is evaluated by including also the localization stage. Fig. 7 shows the scatter plots obtained for the Doppler and Doppler+Doppler-rate based techniques (1000 independent Monte Carlo simulations) for the three cases detailed above: in the sub-figures (b)-(c) the blue/yellow markers pertain the two techniques when using a single antenna receiver, while the orange/violet ones when using the array configuration. For the  same study cases, Table I provides a quantitative analysis showing the standard deviation values of the estimated target velocity components (both methods provide almost unbiased estimates). Starting from ideal conditions [ Fig. 7 (a), case Atarget position exactly known], we can observe as the introduction of the Doppler-rate measure in the estimation procedure yields a significant reduction of the dispersion of the estimated velocity, and hence an improvement of the accuracy. Particularly, Table I, the standard deviation moves from 0.19 m/s and 0.55 m/s to 0.06 m/s and 0.17 m/s for and , respectively. Such improvement is basically maintained also when moving to more realistic conditions, case B in Fig. 7 (b) and case C in Fig. 7 (c), as evident by comparing the yellow to the blue scatter plots and the corresponding values in Table I. Moreover, by comparing the two different receiver configurations, we can also notice that the use of the antenna array receiver [orange/violet markers in Fig. 7 (b) and (c)] allows to recover almost the same performance achievable in ideal conditions [ Fig. 7 (a)]. This is mostly because the target localization is no longer limited to the intersection of the bistatic range ellipses: indeed, the availability of the DoA information largely reduces the localization error thus allowing performance very close to the benchmark even when additional error sources are present.
In order to analyze the improvement provided by the exploitation of the Doppler rate information for different observation geometries, Fig. 8 shows the standard deviation of the estimated velocity components for the same study case of Fig. 7 but by varying the position of the second satellite. Specifically, the performance is analysed as a function of the azimuth separation between the two transmitters being ∆ 1,2 = 2 − 1 and having set 2 = 1 . The same colours coding of Fig. 7 (b-c) is used to characterize the techniques; the markers represent the performance of the two techniques, while the dashed plots the corresponding CRLB values. As it can be noticed:  both the techniques show performance equal to the corresponding CRLB;  for both the velocity components we see a remarkable improvement provided by the Doppler+Doppler-rate based technique, particularly for small values of ∆ 1,2 : such improvement decreases as ∆ 1,2 increases since, in presence of a large angular separation between the bistatic links, the Doppler information becomes suitable for the estimation of both the components;  the best performance is obtained by jointly using the Doppler+Doppler-rate technique and the antenna array configuration. The above analysis proves the effectiveness of the proposed approach in coping with situations involving a not suitable angular diversity between the available transmitters allowing in all conditions the achievement of highly accurate velocity estimates.
The impact of the SNR value on the velocity estimation by means of the Doppler+Doppler-rate based technique is finally examined. As stated previously, the SNR value is taken into account within the weighted estimate (8) by means of the error's covariance matrix Σ . As we already mentioned, in the GNSS framework a large number of satellites characterized by very different values of SNR could be available. In order to test this condition, the previous hypothesis of constant SNR among all the baselines is now eliminated. Specifically, we  and we assume the first two bistatic links sharing the same SNR equal to 1 = 2 = 8 dB and the third one characterized by 3 = 1 ± 10 dB. Fig. 9 shows the performance of the proposed approach with and without including the weighting provided by the error covariance matrix Σ into (8), by considering only the array of antennas configuration. The violet '∇' markers show the results when Σ is included in (8), the yellow 'o' markers when it is omitted and the solid line represents the case when only the two satellites sharing the same SNR are exploited (i.e., for the solid green line = 2). From the figure, we can see that  when the unweighted estimate is used (yellow 'o' markers), the introduction of a third bistatic link causes the performance to improve when SNR3 increases with respect to SNR1 and to degrade with respect to the = 2 case (solid line) when SNR3 decreases with respect to SNR1;  when the weighted estimate is used (violet '∇' markers), as above we observe a performance improvement in the SNR3>SNR1 region and, in contrast, performance is basically maintained equal to the = 2 case when SNR3 decreases with respect to SNR1, thus neglecting the influence of the third bistatic link. This confirms the validity of the proposed approach even in demanding situations where the target response behaves differently on different bistatic links.

VI. EXPERIMENTAL RESULTS
The proposed approach has been tested against real data collected during experimental campaigns carried out in scenarios of practical interest. The experimental hardware was a super-heterodyne receiver collecting signals transmitted by Galileo satellites in the E5a band that transmit pseudo random noise (PRN) codes with chip-rate equal to 10.23 MHz. The reference channel employed a low-gain antenna to record the direct signals from the satellites in visibility, while the radar channel operated with a single element high-gain antenna to collect the weak echoes from the area of interest. The receiving station was also equipped with an Automatic Identification System (AIS) receiver to record in real time the ship actual location, course and speed to be used as ground truth.

A. SCENARIO 1: RIVER SHIPPING
In a first scenario, the passive radar prototype was placed on the bank of the Rhine near Bonn, Germany, observing the river traffic mainly composed by barges and inland cruise ships. Fig. 10 (a) shows the acquisition geometry. During an acquisition of length 60 s, the cargo Addio (49.4 m × 23 m), of which an optical photograph is shown in Fig. 10 (b), was in the field of view of the radar antenna, moving away from the receiver, and GSAT0204 and GSAT0206 satellites, emitting PRN codes E22 and E30, respectively, were acquired. The weakness of the received signals makes hard to detect the target using the conventional MTI processing in Fig. 2. Fig. 11 (a) shows the resulting RD map for GSAT0204 by processing a single frame 3 s long. In the figure, 0 dB denotes the mean background power. As it can be clearly observed, a clear peak corresponding to the target cannot be easily separated from the background, so that the detection is missed with such a short time approach. The long-integration time detection technique in Fig. 3 has been then implemented by considering M = 10 frames of duration 3 s each. Fig. 11 (b) shows the long-time RD map pertaining the Doppler-rate providing the highest gain. It can be seen as the correct TMC allowed concentrating the target energy around the same location, enabling its detection. By applying the long-time integration procedure to both the baselines and shifting each time the integration window of one frame, we achieved 16 data points of the range, Doppler, and Doppler-rate histories for each baseline. As with the particular receiver setup the target DoA was not available, target instantaneous locations were obtained via bistatic range intersection only [5]. Then, the estimated location and the Doppler and Doppler-rate measurements feed the velocity estimation stage.  Fig. 12 shows the estimated velocity components by adopting the conventional Doppler-based method (green 'o' markers) and the Doppler+Doppler-rate approach (red '∇' markers); full line shows the velocity recovered from the AIS message. As it can be observed, in this particular case study the estimates achieved by exploiting the Doppler measurements only are very close to the AIS line. This is mainly due to the fact that the spatial diversity provided by the exploited satellites sufficed to make the linear system (4) wellconditioned. Furthermore, the particular acquisition geometry entailed the targets moving with a dominant radial motion (i.e., low Doppler-rate values). However, from the figures, it is possible to appreciate as, despite the Doppler-rate measurements close to the null value, the proposed approach allows basically maintaining the Doppler-based estimation performance, with a slight improvement. This is also confirmed by the Root Mean Square (RMS) errors of the and estimates reported in Table II, for both methods and taking the AIS information as reference.

B. SCENARIO 2: PORT OPERATIONS
A second acquisition campaign took place near Marghera port, Italy. The receiver was placed at the entrance of the port area, with the radar antenna pointed toward the terminal area acquiring the signals scattered by the commercial vessels entering/exiting the port area, as shown in Fig. 13 (a). Fig.  13 (b) shows an optical photograph of the cargo Tailwind [149.4 m × 23 m], which was moving toward the terminal area during an acquisition of length 250 s. During the acquisition, signals emitted by Galileo satellites GSAT0207 and GSAT0208 (PRN codes E07 and E08, respectively) were collected.
For each baseline, target detection has been performed via the long-integration time approach using the same parameters adopted in the previous case study. By shifting each time the   integration window of 5 s, range, Doppler, and Doppler-rate estimated vectors with 51 elements are obtained.
In the scenario under consideration, ships are constrained to move on a particular direction following the entry/exit route from the port terminal, which is nearly aligned with the radar antenna boresight [see Fig. 13 (a)]. Therefore, the target is also in this case characterized by a dominant radial motion, with Doppler-rate very close to the null value. In order to analyze the impact of the proposed approach against targets characterized by non-negligible values of the Doppler-rate, we artificially rotated the target track by an angle of 45°. To this purpose, the following procedure has been implemented: i. the theoretical range, Doppler, and Doppler-rate tracks are obtained by using the rotated AIS track; ii.
the actual measurements errors are obtained as the difference between the original range, Doppler and Doppler-rate tracks (calculated on the basis of the information recorded by the AIS) and the radar measurements; iii.
finally, the new estimated radar observables are obtained by adding to the theoretical range, Doppler and Doppler-rate tracks, evaluated at step i, the estimation errors pertaining the original track (evaluated at step ii). Fig. 14 (a), (b), and (c) compares the measured target location (obtained by the localization stage), the bistatic Doppler frequency, and the bistatic Doppler-rates, respectively, to the ground truth achieved by the AIS. These measurements are then exploited to perform the target velocity estimation task. Fig. 15 shows the achieved velocity values for the Doppler-based and Doppler+Doppler-rate based estimation techniques. Comparing the results to the AIS ground truth, an improvement of the estimation accuracy can be easily recognized in the first two minutes of the acquisition by exploiting the Doppler-rate information. This is well in line with the theoretical expectations (see Fig. 5), as over such period the experienced Doppler-rate is not negligible, as it can be observed by inspecting Fig. 14 (c). When the Doppler-rate decreases and approaches zero (i.e., the target moves with a dominant radial motion), the two approaches provide similar estimates, still close to the nominal velocity. Table III reports the RMS errors for both the velocity components for the whole  acquisition length (250 s), and the RMS errors of the heading and the velocity magnitude for the first 100 s of the acquisition, showing the improvement provided by the introduction of the Doppler-rate in the estimation procedure.

VII. CONCLUSION
This work introduced a ship velocity estimation technique for passive radar systems relying on satellite opportunistic illuminators and a stationary receiver. The goal is capitalizing on the additional radar measurements made available by the long-integration time procedures usually adopted at the detection stage to overcome the power budget limitations deriving from the choice of using satellite sources. Particularly, an iterative procedure has been implemented exploiting the bistatic Doppler-rates to reach an accuracy gain with respect to conventional approaches relying on the Doppler shifts. The achievable estimation performance has been derived at the analytical level by means of the evaluation of the Cramér-Rao Lower Bound, analyzing the impact of the Doppler-rates exploitation under different system geometries. Few case studies are presented to illustrate the behavior of the method and its effectiveness for different target motion conditions and receivers using a single or an array antenna.
Finally, experimental results using a passive radar prototype and acquiring Galileo satellite signals in different scenarios of interests against targets belonging to different class types have been presented, showing as the approach can be applied with success in practical applications.
While the approach has been here detailed for passive radars based on navigation satellites for maritime surveillance applications, the method could be potentially applied to any multistatic passive radar system, especially in the case of limited power budget conditions. (A10) From (A1) -(A10), it is easy to verify that the FIM can be decomposed as in (11).

ACKNOWLEDGMENT
The experimental campaigns were carried out inside the research project "GALILEO-BASED PASSIVE RADAR SYSTEM FOR MARITIME SURVEILLANCE-SpyGLASS" funded from the European GNSS Agency under