Maritime Moving Target Joint Localization and Kinematic State Estimation Using GNSS-Based Multistatic Radar

Abstract

A global navigation satellite system (GNSS)-based multistatic radar is explored for target localization and kinematic state estimation. Since any point on the earth can be illuminated by a minimum of four satellites of each GNSS constellation at any time, GNSS-based passive radars can be inherently considered multistatic radars. In this paper, a method for jointly estimating the target position and velocity by utilizing both the time delays and Doppler shifts has been proposed, and an analytical accuracy analysis is also provided. In the new method, the bistatic range and Doppler for each path are firstly estimated by the range-Doppler (RD) method, and then by using the bistatic ranges and Doppler shifts. The least squares method is applied to estimate the target position and velocity simultaneously. Compared with the precedent target localization and velocity estimation method, the proposed method achieves a better estimation result with simple procedures. Simulation results are provided to validate the effectiveness of the proposed method.

1. Introduction

Passive bistatic radar (PBR) has attracted a significant growth of attention for target detection in the past two decades. Compared with the active radar, the PBR system has several merits, including covert operations and the fact that it is low-cost and lightweight. Moreover, as it does not emit any electromagnetic wave, the PBR system works silently, which can help avoid attacks via the anti-radiation missile and disturbances via the hostile radar jammer. The PBR system utilizes those illuminators of opportunity, which exist in the environment, so that it does not need the frequency allocations. Numerous research have been conducted for target detection. Among them, a variety of illuminators have been explored, such as FM radio [1,2,3], WIFI [4,5], DVT [6], and so on. The global navigation satellite system (GNSS) [7,8,9,10,11,12,13,14,15,16] is one of the most popular illuminators of opportunity, due to its permanent global coverage, spatial diversity, and plentiful satellite resources (GPS, GLONASS, Galileo, and BeiDou).

This paper focuses on the multistatic radar, based on GNSS signal, for its global coverage and persistent illumination. Any point on the earth can be illuminated by a minimum of four satellites of each GNSS constellation at any time, i.e., a minimum of 16 satellites if all four constellations are considered (GPS, GLONASS, Galileo, and Beidou). It means that the GNSS constellation can be inherently considered as multistatic radar with multiple transmitters and a single receiver. Each satellite and receiver build up a bistatic radar, and the bistatic range and Doppler can be obtained by range-Doppler (RD) processing. Combing the estimated bistatic range and Doppler, the target position and kinematic state can be estimated.

Researchers and engineers have paid much attention to target detection for GNSS-based radar [8,9]. After target detection, estimating the target kinematic state is a natural step. In [7], the author proposed a target localization method, based on the spherical intersection method [17] for GNSS-based multistatic radar. The target position can be well-estimated, and the effectiveness of the method has been verified by real measured data. However, the method only uses time delays and cannot estimate the target’s instantaneous velocity. To estimate the target velocity, a target kinematic state estimation method has been discussed in [14]. The feasibility of the proposed approach, as well as its accuracy, have been discussed, and the target kinematic state has been estimated via experimental results. However, the approach cannot estimate the target’s position and velocity simultaneously. Moreover, the approach needs to obtain the target’s initial position in advance for velocity estimation. Usually, the initial position of the target is unknown, due to its non-cooperation. In this approach, authors set the initial position of the target similar to that estimated by [7]. It can be seen that, within the method, there actually exists a process of error propagation, and the position estimation error will degrade the velocity estimation performance.

To solve this problem, a method for jointly estimating the target position and velocity is proposed in this paper. In the new method, the bistatic range and Doppler for each path are firstly estimated by the RD method. With the estimated bistatic ranges and Doppler shifts, the target position and velocity can be obtained simultaneously by the least squares method. The main difference of the method is that is does not need the target position’s prior information, and the error propagation has been avoided effectively. Meanwhile, the localization and velocity estimation accuracy for GNSS-based multistatic radar are analyzed utilizing the Cramer–Rao lower bound (CRLB) [18]. Moreover, once the estimated values have only bistatic ranges, the new method is still effective, and the target position can be obtained independently. Compared with the precedent target localization and velocity estimation method, the proposed method achieves a better estimation result, with simple procedures. Simulation results are provided to validate the effectiveness of the proposed method.

The rest of the paper is organized as follows. Section 2 derives the target position and velocity estimation process and analyzes its accuracy. Section 3 presents the experimental setup and corresponding RD processing results. The target position and velocity results are presented in Section 4, and Section 5 concludes the paper.

2. Joint Localization and Kinematic State Estimation Using GNSS-Based Multistatic Radar

2.1. System Geometry and Signal Model

Figure 1 is the system geometry of the GNSS-based multistatic radar system with M transmitters and a receiver. O-XYZ defines a Christian coordinate system. Without loss of generality, the receiver is set at the origin of the coordinate system, so the location of the receiver is $RX=\mathit{r}={(0,0,0)}^T$. The receiver has two channels: one is used to record the direct signal from the satellite and then for signal synchronization, the other is applied to record the target echoes from the area of interest. The coordinates of the target and transmitters are denoted as

$$TG={\mathit{x}}_t={(x_t,y_t,z_t)}^T$$

(1)

and

$$T{X}_m={\mathit{x}}_i={(x_m,y_m,z_m)}^T$$

(2)

where m is the subscript of the satellites. Hence, the range from the mth transmitter to receiver can be written as

$$R_{Dm}=‖\mathit{r}-{\mathit{x}}_i‖=\sqrt{x_m^2+y_m^2+z_m^2}$$

(3)

The range from target to mth transmitter can be expressed as

$$R_{Tm}=‖{\mathit{x}}_t-{\mathit{x}}_i‖=\sqrt{{(x_m-x_t)}^2+{(y_m-y_t)}^2+{(z_m-z_t)}^2}$$

(4)

and the distance from target to receiver is

$$R_R=‖{\mathit{x}}_t-\mathit{r}‖=\sqrt{x_t^2+y_t^2+z_t^2}$$

(5)

The bistatic range involved in the path from the m-th transmitter can be written as

$$R_m=R_{Tm}+R_R-R_{Dm}$$

(6)

Assuming that the target is moving with velocity ${\dot{x}}_t={(v_x,v_y,0)}^T$, the Doppler frequency of the target echo of the m-th transmitter is given by

$$f_{dm}=\frac{f_c}{c}({\dot{R}}_{Tm}+{\dot{R}}_R)=\frac{f_c}{c}\left[\frac{{\dot{x}}_t^T({\mathit{x}}_t-{\mathit{x}}_i)}{R_{Tm}}+\frac{{\dot{x}}_t^T({\mathit{x}}_t-r)}{R_R}\right]$$

(7)

where ${\dot{R}}_{Tm}$ is the range rate between the target and m-th transmitter, ${\dot{R}}_R$ denotes the range rate from target and receiver, and $f_c$ is the carrier frequency.

The bistatic range and Doppler shift can be measured by the RD processing for a certain path. In this paper, we suppose that the time delays and Doppler shifts have already been estimated by the RD algorithm. Let $\widehat{R}=({\widehat{R}}_1,{\widehat{R}}_2,\cdots ,{\widehat{R}}_M)$ and ${\widehat{f}}_d={({\widehat{f}}_{d1},{\widehat{f}}_{d2},\dots ,{\widehat{f}}_{dM})}^T$ be the vectors of the estimated bistatic range and Doppler shift.

We now utilize the estimated bistatic range vector $\widehat{R}$ and Doppler shift vector ${\widehat{f}}_d$ to obtain the position and velocity of the target. By rewriting and rearranging Equation (6), we have

$$R_m+{\left[{\mathit{x}}_i^T{\mathit{x}}_i\right]}^{\frac{1}{2}}-{\left[{\mathit{x}}_t^T{\mathit{x}}_t\right]}^{\frac{1}{2}}={\left[{({\mathit{x}}_t-{\mathit{x}}_i)}^T({\mathit{x}}_t-{\mathit{x}}_i)\right]}^{\frac{1}{2}}$$

(8)

After squaring both sides and reformulating them, we obtain

$$R_m^2+2R_m{R}_{Dm}=(2R_m+2R_{Dm}){\left[{\mathit{x}}_t^T{\mathit{x}}_t\right]}^{\frac{1}{2}}-2{\mathit{x}}_i^T{\mathit{x}}_t$$

(9)

Taking the time derivative of (9), we have

$$R_m{\dot{R}}_m+R_m{R}_{Dm}={\dot{R}}_m{\left[{\mathit{x}}_t^T{\mathit{x}}_t\right]}^{\frac{1}{2}}+(R_m+R_{Dm}){\left[{\mathit{x}}_t^T{\mathit{x}}_t\right]}^{-\frac{1}{2}}{\mathit{x}}_t^T{\dot{x}}_t-{\mathit{x}}_i^T{\dot{x}}_t$$

(10)

Equations (9) and (10) are combined to yield the following matrix form:

$$AX=B$$

(11)

where $A$ is

$$A=\left[\begin{array}{cccc}2R_1+2R_{D1}& -2{\mathit{x}}_1^T& 0& {\mathbf{0}}_3^T\\ ⋮& ⋮& ⋮& ⋮\\ 2R_M+2R_{DM}& -2{\mathit{x}}_M^T& 0& {\mathbf{0}}_3^T\\ {\dot{R}}_1& {\mathbf{0}}_3^T& R_1+R_{D1}& -{\mathit{x}}_1^T\\ ⋮& ⋮& ⋮& ⋮\\ {\dot{R}}_M& {\mathbf{0}}_3^T& R_M+R_{DM}& -{\mathit{x}}_M^T\end{array}\right]$$

(12)

$X$ is the parameter related to the position and target velocity of the target

$$X=\left[\begin{array}{c}{\left[x_t^T{x}_t\right]}^{\frac{1}{2}}\\ {\mathit{x}}_t\\ {\left[x_t^T{x}_t\right]}^{-\frac{1}{2}}{\mathit{x}}_t^T{\dot{x}}_t\\ {\dot{x}}_t\end{array}\right]$$

(13)

and $B$ is the constant vector

$$B=\left[\begin{array}{c}{R}_1^2+2R_1{R}_{D1}\\ ⋮\\ R_M^2+2R_M{R}_{DM}\\ R_1{\dot{R}}_1+R_1{R}_{D1}\\ ⋮\\ R_M{\dot{R}}_M+R_M{R}_{DM}\end{array}\right]$$

(14)

where ${\mathbf{0}}_3^T$ is a $3\times 1$ zero column vector.

Then, the target position and velocity can be obtained by calculating (11). The preliminary solution can be expressed as

$$X={(A^{T}A)}^{-1}{A}^{T}B$$

(15)

It is worth mentioning that the proposed method can also estimate the target position when only the time delays are available. This situation can be regarded as a special case, where the Doppler measurements are all zero. Substituting ${\dot{R}}_m=0$ into (12) and (13), the target position can be obtained by only using the estimated time delays.

2.2. Accuracy Analysis

The coefficient matrix $A$ is determined by the satellite positions, estimated bistatic ranges, and Doppler shifts. To analyze the estimation performance of the proposed method, the CRLB is derived. We denote the estimation parameter vector as $\mathit{\theta }={[{\mathit{x}}_t^T,{\dot{x}}_t^T]}^T$ and $\mathit{\alpha }={[R^T,{\dot{R}}^T]}^T$, respectively. The logarithm of the probability density function under $\theta$ can be written as

$$\mathrm{In}p(\mathit{\alpha }|\mathit{\theta })=-\frac{1}{2}{(\mathit{\alpha }-{\mathit{\alpha }}^0)}^T{Q}^{-1}(\mathit{\alpha }-{\mathit{\alpha }}^0)$$

(16)

where ${\mathit{\alpha }}^0$ and $Q$ are the true value and covariance matrix of the measurement vector. The covariance matrix $Q$ can be calculated as

$$Q=\mathrm{diag}({\mathit{Z}}_{M\times M},{\mathit{\Omega }}_{M\times M})$$

(17)

The covariance matrix of the bistatic range and Doppler shifts as

$$\begin{array}{ll}{Z}_{M\times M}& =\mathrm{E}[R{R}^T]-\mathrm{E}[R]{(E[\mathit{R}])}^T\\ & =\mathrm{diag}({\sigma }_{R1}^2,{\sigma }_{R2}^2,\dots \dots ,{\sigma }_{RM}^2)\end{array}$$

(18)

and

$$\begin{array}{ll}{\Omega }_{M\times M}& =\mathrm{E}[\dot{R}{\dot{R}}^T]-\mathrm{E}[\dot{R}]{(E[\dot{R}])}^T\\ & =\mathrm{diag}({\lambda }^2{\sigma }_{fd1}^2,{\lambda }^2{\sigma }_{fd2}^2,\dots \dots ,{\lambda }^2{\sigma }_{fdM}^2)\end{array}$$

(19)

where ${\sigma }_{Rm}^2$ and ${\sigma }_{fdm}^2$, respectively, represent the variance of the bistatic range and Doppler of the mth satellite. It should be noted that the matrix is diagonal. The variance of the bistatic range and Doppler is related to the signal-to-noise ratio (SNR), which can be expressed as [14]

$${\sigma }_{Rm}^2=\frac{c^2}{B^{2}SN{R}_m}$$

(20)

$${\sigma }_{fm}^2=\frac{90}{{\pi }^2{T}_{coh}^2}\frac{1}{SN{R}_m}$$

(21)

By definition, the Fisher information matrix (FIM), whose inverse yields the CRLB, can be calculated as

$$\mathrm{FIM}(\mathit{\theta })=E\left[\frac{\partial \ln p(\alpha |\theta )}{\partial \theta }{\left(\frac{\partial \ln p(\alpha |\theta )}{\partial \theta }\right)}^T\right]={\left(\frac{\partial \alpha }{\partial \theta }\right)}^T{Q}^{-1}\left(\frac{\partial \alpha }{\partial \theta }\right)$$

(22)

where $\frac{\partial \alpha }{\partial \theta }$ is the partial derivative, with respect to target parameters, and can be written as

$$\frac{\partial \alpha }{\partial \theta }=\left[\begin{array}{cc}\frac{\partial R}{\partial x_t}& \frac{\partial R}{\partial {\dot{x}}_t}\\ \frac{\partial \dot{R}}{\partial x_t}& \frac{\partial \dot{R}}{\partial {\dot{x}}_t}\end{array}\right]$$

(23)

The entries of sub matrices are given as follows:

$$\frac{\partial R}{\partial {\mathit{x}}_t}=\frac{{({\mathit{x}}_t-{\mathit{x}}_i)}^T}{R_{Tm}}+\frac{{\mathit{x}}_t^T}{R_R}$$

(24)

$$\frac{\partial \dot{R}}{\partial {\mathit{x}}_t}=\frac{{\dot{x}}_t^T}{R_{Tm}}-\frac{{({\mathit{x}}_t-{\mathit{x}}_i)}^T{\dot{R}}_{Tm}}{R_{Tm}^2}+\frac{{\dot{x}}_t^T}{R_R}-\frac{{\mathit{x}}_t^T{\dot{R}}_R}{R_R^2}$$

(25)

$$\frac{\partial \dot{R}}{\partial {\dot{x}}_t}=\frac{{({\mathit{x}}_t-{\mathit{x}}_i)}^T}{R_{Tm}}+\frac{{\mathit{x}}_t^T}{R_R}$$

(26)

For $m=1,2,\dots ,M$, and $\frac{\partial R}{\partial {\dot{x}}_t}=0$. Then, the CRLB can be given by

$$\mathrm{CRLB}\left(\theta \right)=\mathrm{FIM}{(\theta )}^{-1}$$

(27)

3. Simulation Results

To validate the feasibility and effectiveness of the proposed method, the simulation results were performed with six satellites by using the MATLAB software. The satellite and radar parameters are shown in

Table 1. Simulation parameters of the satellites and signals.

	Parameters	Value
GPS signal parameters	Carrier frequency	1575.42 MHz
	Signal bandwidth	1.023 MHz
	Pulse repetition frequency (PRF)	1000 Hz
Satellite 1	Satellite azimuth angle	37.88°
	Satellite elevation angle	63.44°
	Satellite-to-receiver range	20,387,437 m
Satellite 2	Satellite azimuth angle	112.43°
	Satellite elevation angle	50.10°
	Satellite-to-receiver range	21,638,002 m
Satellite 3	Satellite azimuth angle	214.51°
	Satellite elevation angle	51.65°
	Satellite-to-receiver range	21,302,318 m
Satellite 4	Satellite azimuth angle	312.64°
	Satellite elevation angle	64.21°
	Satellite-to-receiver range	20,657,765 m
Satellite 5	Satellite azimuth angle	233.68°
	Satellite elevation angle	21.53°
	Satellite-to-receiver range	23,620,706 m
Satellite 6	Satellite azimuth angle	291.87°
	Satellite elevation angle	7.58°
	Satellite-to-receiver range	24,981,158 m
Radar Parameters	Sampling frequency	4.092 MHz
	Coherent processing interval (CPI)	1 s

. GPS L1 signal was adopted as the transmitter of opportunity. The initial position of the target was $\left(1000,200,0\right)m$, and the velocity was $\left(6,8,0\right)m/s$. The signal acquisition time was 30s.

To simulate the condition, in line with the experimental conditions, Gaussian white noise was added to the reflected signal as the background disturbance. In the simulation, we only considered the target echo; the direct-path and multi-path signal were not considered. The cancellation of direct-path and multi-path signals has been well-illustrated in the previous work [19,20], and is beyond the scope of this paper.

3.1. CRLB Analysis

Firstly, the CRLB of the proposed method was analyzed, and the signal parameters are shown in Table 1. As seen in Equations (20) and (21), the variance of the bistatic range and Doppler is related to the SNR. The SNR of each target after signal processing is related to the signal processing time and target radar cross section (RCS). As already underlined, different illumination angles, due to the multiple satellites, result in a large variation of the target RCS. Next, the CRLB of the proposed method for two scenarios, considering the different variations of the average RCS observed by multiple bistatic geometries, are provided to analyze the relationship between the SNR, target localization, and velocity estimation accuracy.

Set the corresponding SNR of the target echo of the m-th satellite after signal processing as $SN{R}_m=a_{m}SNR$, where the coefficient $a_m$ meets the requirement of ${\displaystyle {\sum }_m{a}_m}=M$. At first, we consider that the target echo SNR after the signal processing of each satellite is the same, which means that the weight coefficients $a_m=1,m=1,2,\dots ,M$. The localization and velocity estimation performance of the proposed method was analyzed under different noise levels. The SNR $SN{R}_i$ changed from 5 to 25 dB, with a step length set as 1 dB. A total of 100 Monte Carlo simulations were performed under each SNR to obtain the CRLB. Moreover, to analyze the performance of the proposed method under various illumination sources, the CRLB with different satellite number was also simulated.

Figure 2a shows the CRLB of the target localization, with four, five, and six satellites. It can be seen in Figure 2a that, with a higher SNR of the target echo, a better localization performance is guaranteed. Additionally, it was shown that the CRLB of the target position was much lower as the number of satellites increased, which indicates that the target localized more accurately. Once the SNR was low, the variation of the position CRLB for different number satellites was large. Additionally, the variation became very small when the SNR of each satellite increased to more than 20 dB. The CRLB of the target velocity estimation, shown in Figure 2b, had the same estimation performance.

Next, the CRLB of the target localization and velocity estimation, under the condition that the SNR of each satellite is different, was analyzed. Three groups of parameters were simulated, and the parameter settings were $a_{1,2,3,4,5,6}=1$, $a_1=2, a_{2,3}=0.5, a_{4,5,6}=1$, and $a_1=4, a_2=1, a_{3,4,5,6}=0.25$, respectively. Figure 3 illustrates the CRLB of target position and velocity versus SNR under different parameters. Figure 3a shows that, when the SNR of each satellite was different after signal processing, the target localization accuracy was lower than that of equal SNR. Moreover, comparing the CRLB positions of parameter groups 2 and 3, it is seen that a greater SNR variation between each satellite led to a poorer localization performance. The CRLB of the velocity estimation was consistent with the CRLB of the target position.

3.2. Bistatic Ranges and Doppler Extraction

Target localization and velocity estimation are based on the fact that the bistatic ranges and Doppler shifts have been estimated. Next, we will focus on the bistatic ranges and Doppler shifts extraction.

Usually, the bistatic range and Doppler are obtained by the RD processor. At first, a signal synchronization algorithm was applied to track the parameters of the direct signal (delay, Doppler, and phase) from each satellite at the reference channel. With the tracked parameters, a noise-free replica of the direct signal, which was reference signal, was constructed for matched filtering along the range dimension. After that, Fourier transform was applied along the slow time dimension to separate the targets with different Doppler characteristics with appropriate coherent processing interval. Then, a RD map was obtained, and the bistatic range and Doppler of the target could be extract from the RD map. Considering that the target would migrate the range cell and Doppler cell during the coherent processing interval, the Fourier transform, along the slow time dimension in the RD processor, was substituted by the long-time integration Fourier transform (LIFT), which was discussed in detail in [21], so only a brief description is provided here to avoid duplication.

Denote the discrete form of the m-th target echoes as $S_m(p,q)$ and RD map as ${\mathrm{RD}}_m(r,f^i)$, where $p$ and $q$ are denote the fast time and slow time, respectively. Then, the coherent integration process can be expressed as

$${\mathrm{RD}}_m(r,f^m)=F^{\prime }\cdot S_m(p,q)$$

(28)

where $F^{\prime }$ is the LIFT and can be expressed as

$$F^{\prime }=W^{q{q}^{\prime }}\otimes W^{{(q-Q/2)}^2{q}^{\prime }\mu /2/PRF}$$

(29)

where $W=\exp (-\mathrm{j}\frac{2\pi }{Q})$, $q,q^{\prime }=0,1,2,\dots ,Q-1$, $\otimes$ denotes the Hadamard product, $PRF$ is the pulse repetition frequency of the radar system, and $\mu$ is the parameter that is related to the target motion and can be obtained by parameter estimation methods.

The RD map of the target is shown in Figure 4 and consists of the results obtained by six satellites. The RD images indicate that LIFT can well-accumulate the target signal energy, and the target was clearly visible in the RD images. Moreover, it can be seen in Figure 4 that the bistatic range and Doppler for each satellite were different, which was caused by the difference in bistatic geometry. By detecting the peak of the RD map, the bistatic range and Doppler can be obtained. In order to achieve target localization and kinematic state estimation, the bistatic range and Doppler history of each satellite were needed and could be extracted by each RD map at consecutive CPIs.

Figure 5 provides the extracted bistatic ranges versus the time for different satellites. For brevity, the bistatic ranges of six satellites are plotted in a map. It can be seen from the figure that the bistatic distance changed uniformly with the time. Meanwhile, it is seen that, for a certain satellite, different times may correspond to the same bistatic distance, which was mainly caused by the following two reasons: one is that the target is moving slowly, and the other is the poor range resolution. Usually, the range resolution depends on the signal sampling frequency. In the simulation, the signal sampling frequency was set as four times the signal bandwidth. Hence, the range resolution was 73.3 m. In order to improve the accuracy of bistatic distance extraction, we can increase the signal sampling frequency of the radar system.

The bistatic Doppler can also be obtained by the RD map by taking the peak of the target echo. Figure 6 shows the bistatic Doppler variation curve versus time for different satellites. It can be seen from the figure that the Doppler changed uniformly with slow time for the different satellites. Compared with bistatic ranges, the variation of Doppler with time was much more obvious. At the same time, due to the different bistatic structure, the Doppler corresponding to the target under different illumination sources was quite different.

3.3. Localization and Velocity Estimation Results

Figure 7 shows the target positioning results using four, five, and six satellites, respectively. Figure 7a shows the results of target positioning with four satellites. It can be seen from the figure that the localization point was far from the real position of the target, and the localization error was large, up to hundreds of meters at the maximum. When five satellites were used for target positioning, the accuracy of target positioning results was significantly improved, and the positioning points were evenly distributed on both sides of the real motion trajectory of the target. When the number of satellites reached six, the number of localization points distributed on both sides of the target trajectory was more and more dense, which proves that the localization accuracy of the system has been further improved. The localization results of the proposed method were consistent with the localization accuracy analysis of the CRLB. It can be concluded that more satellites can guarantee a higher localization performance, and the localization results by [7] also validate this point.

Figure 8 shows the target velocity estimation results under four, five, and six satellites. It can be seen that when the satellite number was four, due to the initial target position is set as the position estimated by [7], the estimation accuracy in [13] was low. The target estimation accuracy by the new method is gradually improved over the time. When six satellites were considered, the two methods achieved a proximate estimation result. It can be inferred that, with an increase of the number of satellites, the velocity estimation accuracy of the system will approach the real one.

Next, the root means square error $\mathrm{pos}\_\mathrm{RMSE}$ and $\mathrm{vel}\_\mathrm{RMSE}$ were used to quantity evaluate the localization and estimation accuracy of the system. The definition of the $\mathrm{pos}\_\mathrm{RMSE}$ and $\mathrm{vel}\_\mathrm{RMSE}$ can be expressed as

$$\mathrm{pos}\_\mathrm{RMSE}=\sqrt{\frac{1}{I}{\displaystyle \sum _{i=1}^I{({\mathit{x}}_t-{\widehat{x}}_t^i)}^2}}$$

(30)

$$\mathrm{vel}\_\mathrm{RMSE}=\sqrt{\frac{1}{I}{\displaystyle \sum _{i=1}^I{({\dot{x}}_t-{\widehat{\dot{x}}}_t^i)}^2}}$$

(31)

where ${\mathit{x}}_t$ and ${\widehat{x}}_t^i$ represent the real position of the target and the estimated value of the target position, respectively; ${\dot{x}}_t$ and ${\widehat{\dot{x}}}_t^i$, respectively, represent the real velocity of the target and the estimated one, and $I$ represents the number of CPIs.

Figure 9 shows the variation curve of the $\mathrm{pos}\_\mathrm{RMSE}$ with the number of satellites. It can be seen from Figure 9 that, when four satellites were used for target localization, the root mean square error of localization results by [13] could reach 200 m. When the number of positioning satellites reached more than five, the root mean square error of system positioning was within 50 m. Meanwhile, it can be seen that, when the satellite number reached to six, the new method, as well as the method in [7], could achieve the same localization accuracy.

Figure 10 shows the change curve of the $\mathrm{vel}\_\mathrm{RMSE}$ with the number of satellites. It can be seen from the figure that the velocity estimation error of the system was within 1m/s. When the number of satellites of the system reached more than five, the $\mathrm{vel}\_\mathrm{RMSE}$ of the target velocity estimation was within 0.1 m/s, which proves the effectiveness of using the system to extract the target motion parameters.

In the above, the GNSS-based multistatic PBR system was used to locate the target and extract the motion parameters. Simulation results show that the GNSS-based multistatic PBR system can effectively estimate the position and motion information of the target. At the same time, it shows that the accuracy of target positioning and motion parameter extraction of the system were related to the number of satellites. The more the number of satellites, the higher the accuracy of the system.

4. Discussions

Any point on the earth can be illuminated by a minimum of four satellites of each GNSS constellation at any time, and GNSS-based passive radar can be inherently considered multistatic radars. Once the target is detected, position and velocity estimation are a natural step. Joint target localization and kinematic state estimation were proposed in this paper, and now the comparison of the new method with the precedent methods are provided.

Usually, target localization and kinematic state estimation are two separate steps. Considering that the receiver cannot obtain the direction of arrival (DOA) of the target, the position of the target can be obtained with a minimum of four satellites. Moreover, kinematic state estimation always needs the initial position of the target. In the precedent methods, the initial position is set as the estimation position obtained by the localization method. However, the position estimated is usually inaccurate, and there exists a process of error propagation, which will seriously degrade the velocity estimation performance. To improve the accuracy of the target localization and velocity estimation, the new method has combined the separate two together, and the error propagation is eliminated. Hence, better estimation results are guaranteed.

It should be noticed that the new method uses bistatis ranges and Doppler shifts in the estimation process. Once the estimated values have only bistatic ranges, the new method is effective, and the target position can be obtained.

The feasibility and effectiveness of the proposed method has been proven by the simulated data. We achieved the experimental results to validate the effectiveness of the proposed method and explore the capacity of GNSS-based multistatic radar for joint target detection is our future work.

5. Conclusions

A method jointly estimates the target position and kinematic state for a GNSS-based multistatic radar, which has been proposed in this paper. In the new method, the bistatic range and Doppler for each path are firstly estimated by the RD method. With the estimation, as well as by using the bistatic ranges and Doppler shifts, the target position and velocity can be obtained simultaneously by the least squares method, due to the fact that the method does not need the target position prior information, and the error propagation has been avoided effectively. Moreover, the method can obtain the target position independently when the estimated values only have bistatic ranges. Compared with the precedent target localization and velocity estimation method, the proposed method achieves a better estimation result, with simple procedures. Simulation rsults are provided to validate the effectiveness of the proposed method.