Moving Target Detection in Multi-Static GNSS-Based Passive Radar Based on Multi-Bernoulli Filter

: Over the past few years, the global navigation satellite system (GNSS)-based passive radar (GBPR) has attracted more and more attention and has developed very quickly. However, the low power level of GNSS signal limits its application. To enhance the ability of moving target detection, a multi-static GBPR (MsGBPR) system is considered in this paper, and a modified iterated-corrector multi-Bernoulli (ICMB) filter is also proposed. The likelihood ratio model of the MsGBPR with range-Doppler map is first presented. Then, a signal-to-noise ratio (SNR) online estimation method is proposed, which can estimate the fluctuating and unknown map SNR effectively. After that, a modified ICMB filter and its sequential Monte Carlo (SMC) implementation are proposed, which can update all measurements from multi-transmitters in the optimum order (ascending order). Moreover, based on the proposed method, a moving target detecting framework using MsGBPR data is also presented. Finally, performance of the proposed method is demonstrated by numerical simulations and preliminary experimental results, and it is shown that the position and velocity of the moving target can be estimated accurately. in Bayesian recursive and multi-Bernoulli filter.

To overcome the low power level for target detection in a GBPR system, advanced signal processing methods, such as SAR imaging and two-dimensional (2-D) coherent integration [10][11][12][13], are needed to achieve the 2-D coherent integration gain [14,15].Normally, the GBPR system is based on the bi-static model, where a single GNSS satellite is used as the transmitter.However, GNSS constellations (GPS, BeiDou, GLONASS, and Galileo) possess over 100 satellites, and multiple azimuth and elevation angle observations are available for target illumination in any given area [1,2]; on the other hand, the pseudo random noise (PRN) codes used in GNSS have excellent auto-correlation and cross-correlation properties [3].As a result, it is feasible to employ multiple transmitters in a GBPR system at the same time, and the echo signal from each transmitter can be easily separated by the corresponding PRN code.
Our research in this work is focused on moving target detection using a multi-static GBPR (MsGBPR) system, so that information from different observation perspectives can be fully utilized for enhanced detection ability and more effective detection performance.For such a multi-sensor detection problem, centralized, distributed or decentralized architectures are commonly used, and the centralized one is the optimal architecture in terms of detection performance [16].A comprehensive overview of multi-sensor-based detection techniques can be found in [16][17][18] and references therein.The Random Finite Set (RFS) and Finite Set Statistics (FISST) frameworks proposed by Mahler [17,18] have gained much attention and several popular filters such as the probability hypothesis density (PHD) filter [19], the cardinalized PHD (CPHD) filter [20], the multi-Bernoulli (MB) filter [21], the Poisson MB mixture (PMBM) filter [22], and the generalized labeled MB (GLMB) filter [23,24], have been developed based on the FISST framework.These filters also belong to the track-before-detect (TBD) technique [25], which is suitable for low signal-to-noise ratio (SNR) applications, as it minimizes information loss by sidestepping the detection implementation [25][26][27][28][29]. Thus, these filters are good candidates for target detection in MsGBPR with its low power level.The PHD filter and MB filter with subsequent improvements were proposed to cater for multi-sensor applications [30][31][32][33][34][35][36][37][38][39].Compared with the class of PHD filters, the MB-based is more efficient and accurate if each target's existence probability is required [21,39], and the performance of MB-based filters has been well demonstrated in many practical applications, such as computer vision [40], acoustic sensor tracking [41], and sensor control [42,43].
The update strategy is another crucial problem in multi-sensor systems [16,18,22,37,44].In MsGBPR, a single receiver records all target reflected GNSS signals from multiple transmitters, and this can be considered as some kind of centralized architecture.Two well-known strategies for updating in a centralized fusion center are sequential updating and super-measurements updating [18,44].In super-measurement updating, all measurements are first associated to generate a super-measurement, such as the joint likelihood function method; then, the super-measurement will be treated as the measurement by a single sensor [16].In this way, the multi-sensor detection problem is converted into a single sensor detection problem.However, as the joint likelihood function is the product of each sensor's likelihood function, the performance of super-measurement updating is highly dependent on the quality of measurements.In sequential updating, the posterior target state estimated by the measurement from the i-th sensor will be used as the prior information when the measurement from the i + 1-th sensor is being processed [16], such as the iterated-corrected PHD (ICPHD) [18,45] and iterated-corrector MB (ICMB) filters [37].Generally, the sequential updating strategy will yield the best solution, but the estimated result depends on whether the higher quality measurement is used first or last in the sequential updating [18 44].In MsGBPR, the target will be illuminated from multiple azimuth and elevation angles, which will lead to fluctuating SNR levels for different measurements.As the moving target is non-cooperative, the SNR of each measurement is normally unknown.Therefore, it is difficult to decide the update order if the sequential updating strategy is employed.
In this paper, motivated by previous works on 2-D coherent processing, the multi-sensor MB filter is employed in the MsGBPR system for moving target detection.Firstly, the likelihood ratio function in MsGBPR is modeled based on the range-Doppler map.Then, the product MB (PMB) filter and ICMB filter are introduced for target detection in MsGBPR.After that, based on the principle of maximum likelihood estimation, a modified ICMB filter with an online SNR estimation and its sequential Monte Carlo (SMC) implementation are proposed, which can update all measurements in the optimum order and get the best output from the ICMB filter for detection.Finally, the whole moving target detection framework in MsGBPR is presented, including echo separation, 2-D coherent integration, target state estimation by the proposed filter, and target track generation.
This paper is organized as follows.Geometry of the MsGBPR system and its 2-D coherent integration result in range-Doppler domain are presented in Section 2.Then, in Section 3, the modified ICMB filter using MsGBPR data is derived in detail, with its sequential Monte Carlo (SMC) implementation.Furthermore, the moving target detection framework for MsGBPR system is also introduced.In Section 4, numerical simulations are carried out to demonstrate the validity and feasibility of the proposed modified ICMB filter.Discussion about the performance and the computational efficiency improvement by the modified ICMB filter are given in Section 5.In Section 6, preliminary experimental results using GPS signals are also presented to demonstrate the effectiveness of the proposed modified ICMB filter.Section 7 concludes the paper with possible future research directions.

MsGBPR System and Its Signal Model
First of all, the general geometric configuration of an MsGBPR system is described and the echo model for a moving target is derived based on the GNSS signal and MsGBPR acquisition geometry.Then, the map results of a moving target in the MsGBPR system are presented based on 2-D coherent integration processing.

Geometry of MsGBPR and Moving Target Echo Model
Figure 1 shows the general geometric configuration of an MsGBPR system for moving target detection.For MsGBPR, multi-angle observations from over 20 satellites (GPS, BeiDou, GLONASS, and Galileo) in any given area are attainable, and to simplify, only one receiver is employed, which is located at the origin of the Cartesian coordinate system, as shown in Figure 1, where the X-axis points to the North direction and the Y-axis points to the West direction.The moving target P is located at Ptar(XA, YA, ZA) with a constant velocity Vtar(Vx, Vy, 0).Furthermore, multiple global navigation satellites are used as transmitters denoted by T i , as shown in Figure 1, where i = 1, 2, 3, ….For each transmitter T i , its position vector and velocity vector are represented by i Tran

P and i Tran
V , respectively, and both can be obtained precisely using the ephemerides of GNSS satellites.

Moving Target
Receiver The radar system's range history is a crucial factor for signal processing and target detection.Considering the link of the i-th transmitter, target and receiver, the range history Ri(も) can be modeled as [11,12] () where ゆ is the signal wavelength, も is the illumination time, Rref,i, fd,i and fr,i denote the reference range, Doppler centroid, and Doppler modulated rate with the geometry of the i-th transmitter at the center illumination time, and their detailed calculation method can be found in [14].Then, according to [14], after quadrature demodulation and synchronization, the 2-D form of the echo signal can be written as where t is the slow time, わ is the fast time, AT is the signal amplitude, Qi(•) is the modulation code of the i-th transmitter, c is the speed of light, and めi is the target radar cross-section (RCS) associated with the geometry of the i-th transmitter.It should be noted that ゎi is changing with the observation angle, so it may be different for each transmitter geometry.As in (2), the received echo signal of MsGBPR is actually the summation of all reflected signals from visible transmitters (satellites), and it is easy to separate them by each satellite's PRN code, as it has excellent auto-correlation and cross-correlation properties.Unlike the "stop-and-go" approximation in traditional radar [46], the moving of transmitter during the fast time in GBPR should also be considered.Hence, as shown in (2), the Doppler effect during the fast time needs to be taken into consideration, and it should be carefully dealt with during the 2-D coherent integration processing [14].

2-D Coherent Integration Map Results of Moving Target
As the power level of MsGBPR echo signals is very low, long time 2-D coherent integration processing is necessary to achieve sufficiently high SNR for effective target detection.Based on the principle of Radon Fourier Transform (RFT), a modified RFT method is proposed in [14], with range-walk removal, RFT, range compression and phase error compensation.A fast implementation with Chirp-Z transform is also presented in [14].Based on the modified RFT, the diagrammatic sketch of 2-D coherent integration processing and imaging results by MsGBPR are presented in Figure 2. As shown, all visible reflected echo signals are recorded, and multiple range compression results can be obtained using the PRN code of each satellite.Then, moving targets will be focused in the range-Doppler domain by azimuth coherent processing.The range history and Doppler frequency of each geometry configuration (different transmitter used) are totally different.Therefore, as shown in Figure 2, the same target can be located in different range and Doppler cells, which will increase the difficulty of multi-target detection.According to [14], the moving target in GNSS-based radar can be focused in range-Doppler domain, which is ,e x p (3) where QC is the gain of range compression, D(•) is the normalized envelope of the auto-correlation result of the i-th transmitted signal, D(0) = 1, i dv f , is the Doppler frequency caused by target motion, Tp is the signal period, and NA denotes the azimuth sample number.It should be noted that the reference range Rref,i, and target Doppler frequency fdv,i will change as a different geometry with the i-th transmitter is used.
Clearly, the maximum value will occur at the position of target motion parameter (Rref,i, fdv,i), and as in (3), ( ) , , which equals the coherent integration gain.With the method in [14], the 2-D coherent integration gain of several even tens of seconds of echo data can still be achieved.

Moving Target Detection Based on Modified Multi-Bernoulli Filter
In this paper, the MsGBPR system consists of one receiver and multiple transmitters, and the echo signal from each transmitter can be separated by the PRN code.Therefore, the centralized architecture is easy to implement in MsGBPR, and the receiver is the central fusion center.As multiple transmitters are used in MsGBPR, as many as I (number of observable satellites) maps of moving targets can be acquired at the same time, which can be used to improve the target detection performance of the system.
To fully utilize all measurements (map results), multi-sensor target tracking solutions are needed.The multi-Bernoulli filter is a famous tool for multi-sensor systems, which is based on the finite set statistics framework, and it has been widely researched in applications like maneuvering targets, field robotics, and computer vision.Furthermore, the multi-Bernoulli filter is a kind of TBD method, it sidesteps the detection by leveraging spatiotemporal information directly from the map results, and is capable of detecting a moving target with SNR as low as 7 dB [29].Therefore, less coherent integration time is needed or its detection range will increase significantly if the multi-Bernoulli filter is employed.
In this section, both the PMB filter and the ICMB filter for MsGBPR target detection are introduced.To update the multiple map results in the optimal order, a modified ICMB filter with online SNR estimation is proposed.Based on this, a framework for moving target detection in MsGBPR is also presented.

Multi-Bernoulli RFS and Multi-Bernoulli Filter
A multi-Bernoulli RFS X is a union of a certain number of independent Bernoulli RFSs, which can be expressed as = XX  .Its existence probability is r q 樺 (0, 1), and its probability density is p q (•).The probability generation functional of a multi-Bernoulli RFS is [21,37] ( ) . Therefore, a multi-Bernoulli RFS can be represented by the parameter set =  , and the probability density is [21,36] () ( ) Based on the Bayesian recursive framework, the process of a traditional multi-Bernoulli filter can be described as follows.
Prediction: The posterior target density is a multi-Bernoulli parameter set Then, the predicted target density is also a multi-Bernoulli RFS, which is [27] ( ) where is the single target transition density at time k, ps,k is the target existence probability at time k, and the second part in ( 6) is the multi-Bernoulli RFS of newborn targets at time k.
Update: The predicted target density is a multi-Bernoulli parameter set Then, the posterior target density can be approximated by a multi-Bernoulli RFS [27] ( ) where Xk is target state RFS, Zk is the measurement RFS, and L(•) is the likelihood ratio function.
The relationship between Bayesian recursive and multi-Bernoulli filter and their propagation track table are shown in Figure 3.As shown, the target posterior density and target Bayes recursion are performed in a Bayesian filter, and based on RFS, a multi-Bernoulli parameter set is used in prediction and update.

Multi-Beroulli Filtering for Target Detection in MsGBPR
Based on the introduction in Section 3.1, the multi-Bernoulli filter application for target detection in MsGBPR is carried out.Firstly, based on the range-Doppler imaging results, the target model, measurement model and likelihood ratio function are provided.Then, the PMB filter and ICMB filter are presented based on the likelihood ratio function.

Target Model, Measurement Model and Likelihood Ratio Function in MsGBPR
In this paper, targets are assumed to move at a constant speed in the plane, and the target state vector is given by x = [A, x, y, vx, vy] T , where A denotes the amplitude, x and y denote the target positions, and vx and vy are its velocities.Therefore, the traditional target kinematic model can be used for the description of target state transition.
Suppose there are N(k) targets at time k, and their target states are xk,1, …, xk,N(k).Then, the target states RFS is  .Each target xk,qbelongs toXk at time k, either continues to exist at time k + 1 with probability ps,k,q(xk,q) and new state xk+1,q, or dies with probability 1ps,k,q(xk,q).Consequently, for a given target xk,q at time k, its behavior at time k + 1 can be modeled by the Bernoulli RFS Sk+1|k(xk,q).If this target still exists, Sk+1|k(xk,q) is { xk,q }, and if this target dies, Sk+1|k(xk,q) becomes empty set.In addition, new targets could appear at time k + 1.Therefore, target state Xk+1 consists of two parts: existing targets since the last time and new emerging targets, and the RFS Xk+1 at time k + 1 is given by the union [27] ( ) where For measurements or observations in a MsGBPR system, the moving target is focused in the range-Doppler plane, and multiple imaging results can be obtained at each time.Normally, each moving target has a different Doppler frequency and range in MsGBPR, which means there is only one moving target in a pixel unit of the range-Doppler imaging result with no overlapping with each other.Then, considering whether the target exists or not at time k, the map measurement of MsGBPR with the i-th transmitter is given by where Zk,i is the i-th measurement at time k, vk,i is the complex Gaussian noise, with mean value 0 and variance where Ak,i = ATめk,iQC, Rref,i = ぼi(xk), fdv,i = ぽi(xk) and detailed calculation can be found in [14].Based on the above models, it is reasonable to assume that the intensity of each unit follows the Rice distribution if target exists, or Rayleigh distribution if noise only [29,47].Therefore, the likelihood function for target existing and target dying in unit (m, n) is, respectively  (14) where Nf and Nr denote the size of map result in Doppler and range dimensions, respectively.Therefore, for the i-th transmitter, the corresponding statistical likelihood ratio is ( )

PMB Filter
As moving target detection in MsGBPR is a multi-sensor target tracking problem, it is crucial to find a solution to use all observations (maps based on different transmitters).As the transmitters are conditionally independent, the joint likelihood ratio function can be modeled by a product [36] ( ) ( ) Notice that ( 16) has the same functional form as the single-transmitter likelihood ratio function in (15).Therefore, based on the principle of the multi-Bernoulli filter, the joint product likelihood ratio function L(Zk|Xk, ek) can be used for updating, and this kind of multi-sensor multi-Bernoulli filter is called the PMB filter.
In Figure 4, the propagation track table of PMB filter is presented.As shown in Figure 4, the joint likelihood ratio function is a product over all of the likelihood functions, and all measurements at time k are used at the same time, which means that only a single update is needed at time k.Thus, the structure of PMB filter for multiple sensors is the same as the multi-Bernoulli filter for a single sensor, and only a little extra computation load is needed in the PMB filter, when the joint likelihood ratio function is calculated.

ICMB Filter
The PMB filter provides a good solution to the multi-sensor tracking problem, and its structure appears to be potentially computationally tractable.However, if the i-th map result is really poor (SNR is low), the likelihood function Li(zk,i|xk, ek) will be much smaller than the others.In this situation, the joint likelihood function L(Zk|Xk, ek) is mainly decided by this poor quality measurement, which will affect the accuracy of target detection and tracking.
Therefore, the ICMB filter is presented by updating the measurement sequentially.It does not possess a rigorous mathematical framework but provides a practical implementation for multi-sensor tracking [37].As shown in Figure 5, all measurements are used through iterative single map updates.Unlike the PMB filter, a single prediction and multiple updates are required at time k in the ICMB filter.At time k, the posterior parameters are computed by updating the prior parameters with the 1-st map, then treating the updated parameters as the prior parameters and updating those parameters with the 2-nd map, and so forth.This process is repeated until all map measurements are exhausted.Therefore, if the MsGBPR system consists of a large number of transmitters, the computation load will increase significantly.

Modified ICMB Filter with SNR Online Estimation
According to [44], the order of updates plays a major role in the ICMB filter, if the poor quality results and good quality results are both present, and the best strategy is to update the poor quality results first and the good quality results later.However, it is difficult to obtain map quality information before target detection in an MsGBPR system, as the RCS of a target varies significantly with the bi-static angle or observation angle.
Therefore, to get the best target detection result in MsGBPR, the map SNR should be estimated first based on the measurement, and then we can update results iteratively in the order of SNR from the lowest to the highest.The online SNR estimation method is proposed in Section 3.3.1,and then the SMC implementation of modified ICMB filter is presented.

Online SNR Estimation
As shown in (10), the map SNR is decided by the value of amplitude Ai, as the power of noise is already known.Ai varies with the transmitter.To estimate the amplitude information, a traditional way is to increase the target state dimension, such as xk = [xk, yk, vx,k, vy,k, Ak,1, Ak,2, ... , Ak,I] T or xk = [xk, yk, vx,k, vy,k, ゎk,1, ゎk,2, ..., ゎk,I] T .However, this will increase the computation load accordingly, which is not desirable.Moreover, compared with the amplitude information, we are more concerned about the position and velocity of the target, and the amplitude information is only needed when the likelihood function is calculated in the multi-Bernoulli filter.Therefore, it will be much more efficient if the target state is modeled as xk = [xk, yk, vx,k, vy,k] T .
An approximation expression of the Bessel function is [48] () Based on this, the likelihood ratio function with the i-th transmitter can be rewritten as  (18) where . Based on the principle of maximum likelihood estimation, the logarithm is applied on the both sides of ( 18), and the partial derivative of Ak,i is Now, using the estimated ()

SNR
x , we can update the map results in the optimum order in the ICMB filter, and obtain the optimum output for target detection.

SMC Implementation of Modified Iterated-Corrector Multi-Bernoulli Filter
In the following, a modified ICMB filter is proposed based on the online SNR estimation, and a general SMC implementation of the proposed filter is also presented, as shown in Figure 6, where five stages are included: prediction, online SNR estimation, updating, pruning and merging, and target state estimation.For details of SMC implementation of the traditional multi-Bernoulli filter, the reader is referred to [27].

 Prediction
The first stage is prediction, which includes two parts: new birth targets and existing targets.
Assuming the density of birth targets is ( ) , where q ,k r Γ is given by the birth model, and the density q ,k p Γ can be described by the birth particles set , where is the number of birth targets, () , q k L Γ is the number of particles for each birth target, and being the proposed density of birth targets.Normally, uniform distribution is used for the proposed density in the map area, and in this case, the weight of each particle is the same, which is Assuming the predicted density is ( ) , and both possibility ,| 1 q pkk r − and density ,| 1 q pkk p − can be described by the existing particles set where the predict particle state () , ,| 1 jq p kk− x is generated by the target kinematic model using the particle state at time k-1, and the weight of the particle is () ()() j,q j,q j,q p,k|k -1 p,k -1 p,k -1 s,k w= r w p (24) Finally, the predicted targets density, as shown in (6), and the fused particles set , x can be obtained by the particle state, existence probability.and particle weight of the birth particles and existing particles.

 Online SNR estimation
As described in Section 3.2.3, the map SNR can be estimated by (20), and here the maximum likelihood estimation is performed based on the particles set , . For each particle, the amplitude is estimated by (19), and an estimated amplitudes set is obtained based on all particles.It should be noted that, for each map, an estimated amplitudes set can be obtained, with i = 1, 2, ..., I.Then, the estimated amplitude of the i-th map is given by , and the corresponding SNR can be also estimated online, which is Based on ( 25), the map SNR vector is Then, the elements of the SNR vector are arranged in ascending order, and a new vector SNRk is obtained.The ascending order SNRk can be used for the updating stage, and the optimum estimation output is then obtained.

 Multiple updating
Based on (7), the updating stage is performed and the updated posterior targets density Assuming the map SNR in ascending order is [1, 2, …, i,…, I], we need I updates at time k.Then, based on the likelihood ratio function in (14), the i-th updated density is computed as wrw r= p= 1-r +r w w (26) where w= wL | , e r = r w= w zx (27) As shown in ( 26) and ( 27), the (i-1)th updated parameters are treated as the prior parameters for updating those parameters such as ,1 − q ki r and () ,1 − j,q ki w with the i-th map.It should be noted that for the first update (i = 1), the predicted parameters such as q k|k-1 r and () from the prediction stage are used in (26).For the last update (i = I), the parameters ( )  5, after all maps are updated, the modified filter will move to the next stage.

 Pruning, merging, and resampling
Similarly to the traditional multi-Bernoulli filter, pruning, merging, and resampling will be performed after updating.Detail of those implementations can be found in [17,27].The pruning and merging implementation is to reduce the growing number of particles and target tracks, and targets with existence probabilities below the threshold are discarded.To alleviate the degeneracy issue, a resampling process is always performed.For each hypothesized target, particles are resampled, and the number of particles is reallocated.The reallocated number is in proportion to the probability of existence q k r , but it should be larger than the minimum number Lmin per hypothesized track.Thus, after resampling, the reallocated number for the q-th hypothesized target is ( ) where Lmax is the pre-set maximum number of particles per hypothesized track.Finally, the particles set after pruning, merging, and resampling is obtained as

 Target state estimation
The last stage is target state estimation.Normally, the estimated number of targets is . Then, reorder the existence probabilities in descending order to obtain the corresponding position y(idx) of the first ˆk N existence probabilities.Finally, the estimated target states are .More detail about target state estimation can also be found in [18,27].

Moving Target Detection Framework
As shown in Figure 7, the moving target detection framework is presented in this section, which includes signal separation and range compression, 2D coherent processing, target detection based on the multi-Bernoulli filter, target state estimation, and target track estimation.For an MsGBPR system, the receiver can record the echo signal continuously, as shown at the left part of Figure 7. Based on the relationship of coherent time and power budget, the slow time or coherent time T can be decided, which is also the time interval of the map.It should be noted that the time k in Sections 3.2 and 3.3 actually means the slow time from (k-1)T to kT.Then, for each echo signal, signal separation and range compression is performed based on the transmitter's PRN code.As a result, the reflected echo signal from a specific transmitter can be separated and the corresponding range compression result can be also obtained.This is followed by the 2D coherent processing based on modified RFT, and multiple range-Doppler map results are then obtained.After that, depending on the quality or SNR of obtained map results, if the SNR of each map is similar, the traditional PMB filter with likelihood ratio function in ( 16) is implemented for target detection; on the other hand, if the map SNRs are different or unknown, the proposed modified ICMB with online SNR estimation should be adopted, and we update map results iteratively in the ascending order of map SNR.In practice, the modified ICMB filter is recommended, as the map SNR with non-cooperative targets is usually unknown.Then, the updated target density and estimated target states will be the input for the next time instant.Finally, based on the estimated target states at different time instants, targets will be detected and their tracks are also acquired.

Experiments and Results
To demonstrate the performance of the proposed filter and the feasibility of target detection in an MsGBPR system, numerical simulation results including online SNR estimation are provided.Signals from three GPS satellites (Space Vehicle Number, SVN #2, #10, #18) are employed.Assuming the target moves on the XOY plane, the target motion parameters and the associated simulated parameters are listed in Table 1.As listed, 20 maps with each transmitter are simulated, and the coherent time of each map is 2 s.The target appears at frame-5 and disappears at frame-16 with a constant velocity V (-30, 100, 0) m/s.Here the multi-transmitter index is i = 1, 2, 3, and the time series index is k = 1, 2, …, 20.The target motion is modeled by a linear stochastic process where F is the state transition matrix, With T being the frame period, and in our simulations, T = 2 s. vk is the Gaussian distributed noise with zero mean and covariance matrix Q [29].Map formation method in [14] is used for generating the range-Doppler map results, and the map signal model is shown in (11).The SNRk,i is set to 8 dB, 9 dB, 10 dB, with i = 1, 2, 3 corresponding to SVN #2, #10 and #18, respectively.In order to show the range-Doppler maps with different geometric configuration, map results at frame-5, frame-10 and frame-15 using the transmitter signals SVN #2 and #10 are shown in Figure 8.It should be noted that the direct signal in the range-Doppler map has been removed [49,50].To highlight the target, all maps in Figure 8 are noise-free.As shown in Figure 8a,d, the range history of the moving target is totally different and the target will be located in a different range-Doppler cell, if different transmitters are used.Therefore, it is really hard to fuse the range-Doppler map results directly.
To testify the effectiveness of the proposed filter for multiple transmitter signals, target detection experiments based on the set of simulated range-Doppler maps are carried out.For the SMC implementation of the modified ICMB, the birth process is a multi-Bernoulli density The estimated SNR results based on the proposed online estimation method are presented in Figure 9. Based on the simulated map set above, the ideal SNR (SNRk,i) with each transmitter is 8, 9, 10 dB, as shown by the dotted lines in Figure 9. Initially, there is no target from frame-1 to frame-4, so that the estimated SNR are just random values around 0 dB.Then, the target appears at frame-5, and the SNR of map-i is estimated.However, due to the approximation in (18), the estimation accuracy of 5, ML i SNR is a little bit low at the beginning.After that, with more map frames available, the estimation accuracy of , ML ki SNR is gradually enhanced, especially in the high map quality case, which can be seen from the green line in Figure 9. Hence, the performance of the proposed online SNR estimation is better when the map SNR is higher.Moreover, due to the accumulated error of the likelihood ratio function among those map frames, the estimated SNR may be better than the simulated SNR in some frames.However, as shown in Figure 9, the overestimation of map SNR can be overcome as more map frames are used.From frame-16, the target disappears and the estimated SNR becomes randomly distributed again.On the whole, although the estimation accuracy of map SNR is not high at the beginning, the map SNR can be estimated effectively and it will get close to the ideal result, with more and more map frames.Based on the estimated map SNR, the optimum update order (ascending order) is then obtained and the true target track and the estimated one in 2-D space are shown in Figure 10.At the beginning, the map SNR cannot be estimated correctly, causing errors in the modeling of the likelihood ratio function in (15).Thus, the estimated target position is far away from the true position.Then, with more frames, the estimated track becomes close to the true one gradually, demonstrating clearly the effectiveness of the proposed filter.

Discussion
In this section, the performance of the proposed modified ICMB filter is evaluated firstly, in terms of detection possibility, root mean square (RMS) position estimation error, and RMS velocity estimation error.Then, improvement on the efficiency of the proposed filter is presented, with tradeoff between detection efficiency and accuracy.

Performance Evaluation
With the parameters listed in Section 4, firstly, to demonstrate the importance of update strategy, the ICMB filter with both ascending and descending update strategies is used in the experiments.Here, we assume that all map SNRs are already known, so the estimation error in map SNR is avoided; for the purpose of comparison, all map SNRs are either 10 dB or 8 dB.The detected frame rate is defined to represent the detection possibility in Monte Carlo trials, which is Nd,k/Nmc, where Nd,k denotes the number of detected frame in frame-k, Nmc denotes the number of Monte Carlo trials.In our experiments, 1000 Monte Carlo trials are carried out (Nmc = 1000), and the detected frame rate results are shown in Figure 11, with a threshold of 0.5 to determine whether the target is detected or not.As shown, the target is detected when it appears (frame-5), irrespective of the update order.However, the detected frame rate in the ascending update (blue line) case is much higher than that in the descending case (red line), which verifies that the ascending update strategy should be adopted in iterated-corrector filter for multi-sensor detection.Generally, in the iterated multi-sensor filter, the overall performance of the filter will be degraded if the detection probability of the last sensor in the iterated update is low.In radar systems, the detection probability is directly related to the map quality or map SNR.Therefore, the ascending update strategy is the best in the ICMB filter.Moreover, as shown in Figure 11, the detected frame rate in both cases is higher than the case of all map SNRs being 8 dB (green line), but lower than the case of all map SNRs being 10 dB (pink line).Obviously, as expected the high SNR map provides more useful information for target detection.
To show the importance of update order in ICMB filter further, the RMS estimation error of target position and velocity is presented in Figure 12.As shown, both errors are large, but they become smaller with the increase of map frames.Same as the detected frame rate, the RMS error with ascending update order is much smaller than that with descending order, and the RMS errors for both cases are also between the two extreme cases with all-8 dB and all-10 dB map SNRs.Furthermore, as shown by the red line in Figure 12a, the final RMS error of target position in ascending update order is close to 150 m, which is an acceptable RMS position error as the range resolution of MsGBPR is expected to be around that value.Next, the detected frame rate and RMS errors of the ICMB filter with ascending update order and the PMB filter with all map SNRs being modeled as 8 dB and 10 dB are presented.As shown in Figure 13, a missing alarm will appear if all map SNRs are modeled as 8 dB, because the map quality is underestimated; on the other hand, a false alarm also appears when all map SNRs are modeled as 10 dB, as the map quality is overestimated.Obviously, the detected frame rate and RMS error of the ICMB filter and the proposed filter are between the two extreme cases of all 8 dB and all 10 dB, as shown in Figures 13 and 14.However, when the target appears initially, the performance of the ICMB filter is much better than the proposed filter r, as shown by the blue line and red line in Figures 13 and 14.The reason is that the estimation accuracy of map SNR by the proposed estimation method is not high enough, which affects the order of updates and the likelihood function in (15).As shown in Figure 9, the estimation accuracy of map SNR is gradually enhanced with the increasing number of map frames.After a few frames, the map SNR is estimated more accurately, and the ascending update can be realized and the likelihood function can also be modeled accurately.Therefore, as shown by the blue line and red line in Figures 13 and 14, the performance of those two filters is almost the same after a few frames, and the effectiveness and accuracy of the proposed filter are also verified.

Improvement in Computational Efficiency
For the multi-Bernoulli filter, one of the time-consuming processes is the calculation of existence probability q k r in (8), which depends mainly on the calculation of the likelihood ratio in (15).In the MsGBPR system, the transmitted code in each GNSS satellite is the high-rate PRN sequence, which has a particularly excellent auto-correlation property.So, for the resultant range-Doppler map, the peak value and energy of side-lobe in range direction are sufficiently small.In addition, the azimuth envelope of the map is in the form of a sinc function.Therefore, the target energy is mainly concentrated in the main lobe, and the sub-area can be substituted for the whole area to calculate the likelihood ratio, which can be rewritten as  (30) where Nf_s and Nr_s denote the size of the sub-area map in Doppler and range domains, respectively.
Experiments with different sizes of sub-area are carried out in this section.Seven different sizes are adopted for the calculation of the likelihood ratio, which are 1 × 1, 3 × 3, 5 × 5, 7 × 7, 9 × 9, 11 × 11 and the whole map size, and the corresponding performance metrics of the proposed filter are listed in Table 2.The power ratio in Table 2 denotes the energy of the sub-area divided by the energy of the whole area.As shown, almost 90% energy is concentrated in the 5 × 5 sub-area.The corresponding computation time and detected frame rate result for different sizes of sub-area are also compared in Table 2.As anticipated, the substitution of the sub-area for the whole map to calculate the likelihood ratio will improve the computational efficiency significantly.The smaller the sub-area size is, the less the required computation time.However, substituting a single cell for the whole map to calculate the likelihood ratio improves the efficiency most, but it degrades the detected frame rate seriously.In addition, the performance degradation of detected frame rate also takes place in case 2 (3 × 3), but is negligible for other cases.Over all, the 5 × 5 sub-area is a better choice as a substitution for the whole area, with a good tradeoff between accuracy and efficiency.With the recorded moving target echo signal, 11 frames of range-Doppler map are generated for each transmitter (SVN #12 or #25), and the time interval or coherent integration time is 0.1 s.Since the coherent integration time is very short, the airplane target can be assumed to move at a constant speed.Moreover, as the average vertical speed of airplane in a normal landing is around 2 m/s, the airplane can also be assumed to move in the XOY plane.However, as shown in Figure 16, the SNR of range-Doppler map is very low, and it will be difficult to tell whether the moving target exists directly from the range-Doppler map.In Figure 16, both the Doppler frequency caused by the satellite moving and the reference range are removed.Then, based on the proposed modified ICMB filter, the moving target detection experiment is carried out.According to the experimental scene in Figure 15, the target mainly moves in the X direction, and normally, the speed of a landing airplane is around 75 m/s.So, the initial state of the target can also be generated  17, and the map SNR using SVN #12 and #25 is around 11 dB and 12 dB, respectively.As demonstrated in Figure 17, the map SNR can be estimated accurately by the proposed online estimation method.So, the map using SVN #12 should be updated first during the sequential updating.In the detection experiment, the moving target is detected in 9 map frames, and it moves from frame-2 to frame-10 with the estimated speed vector being (72.4,3.2) m/s.The main reason that the target is not detected in frame-1 and frame-11 may be due to the target being out of the receiver antenna beam.The estimated track is presented in Figure 18.The estimated speed scalar is 72.5 m/s and the estimated target range is about 200 m, which are reasonable results demonstrating again the effectiveness of the proposed modified ICMB filter.

Conclusions
In this paper, an MsGBPR system with a single receiver has been studied for target detection, and a modified MB filter with online SNR estimation employed as the MB filter is well suitable for low SNR applications and has good performance in multi-sensor tracking.Since the RCS of a target varies significantly with the bi-static angle or observation angle, the resultant map SNRs for different transmitters are different and both sequential updating and super-measurements updating strategies are difficult to implement in a traditional MB filter.In the proposed solution, first, the likelihood ratio function is modeled based on the range-Doppler map; based on the maximum likelihood estimation method, an online map SNR estimation method is developed using the map measurements and target signal model; then, the modified ICMB filter is proposed, which can update all measurements for multi-transmitters in the optimum order (ascending order).Furthermore, a moving target detecting framework for MsGBPR is also presented, including echo separation, 2-D coherent integration, target state estimation by the proposed method and its track generation.The excellent performance of the proposed method has been demonstrated by numerical simulations and some preliminary experimental results.Moreover, to improve the computational efficiency, a sub-area can be used as a substitute for the whole map to calculate the likelihood ratio, and based on simulation results, a 5 × 5 sub-area is recommended for the proposed method, which provides a good tradeoff between computational efficiency and detection accuracy.
Based on findings in this paper, our next work is to test the proposed method for multiple moving targets with real MsGBPR data, and later work will explore target detection using multiple

Figure 1 .
Figure 1.General geometric configuration of a multi-static global navigation satellite system based passive radar (MsGBPR) system.

Figure 2 .
Figure 2. Diagrammatic sketch of 2-D coherent integration processing and map results in MsGBPR.

Figure 3 .
Figure 3. Propagation track table in Bayesian recursive and multi-Bernoulli filter.

Figure 4 .
Figure 4. Propagation track table of the product multi-Bernoulli filter.

Figure 6 .
Figure 6.Sequential Monte Carlo (SMC) implementation of the proposed modified ICMB filter.
. Unlike the traditional multi-Bernoulli filter, the iterative updates are implemented based on the estimated SNRk obtained from the online SNR estimation stage.
estimated target RFS is { } ˆ1 ˆk N y k y= x

Figure 7 .
Figure 7.The proposed moving target detection framework for the MsGBPR system.

Figure 10 .
Figure 10.Target detection results by the modified ICMB filter.

Figure 11 .Figure 12 .
Figure 11.Detected frame rate in different update order cases.
by the uniform distribution, where x~U[ 200, 200 m], y~U[50 m, 300 m], vx~U[60 m/s, 80 m/s], vy ~U[0 m/s, 20 m/s].In addition, those parameters of the multi-Bernoulli filter are the same as those listed in Section 4. Using the online SNR estimation method, the estimated SNRs and measured SNRs are shown in Figure

Figure 17 .
Figure 17.The estimated SNR results using the experiment data.

Figure 18 .
Figure 18.The estimated target track using the experiment data.
(3)k is a target model variable indicating whether the target exists at time k or not, and gk,i(xk) is the target response function, as presented in(3).Moreover, target motion parameters are used in the response function, but the target state vector is x = [A, x, y, vx, vy,] T .As a result, the target reference range and Doppler frequency should be calculated first, based on the predicted or estimated target state x, and the response function gk,i(xk) can be rewritten as is the zero-order Bessel function.Based on the ambiguity function of GNSS signals, the correlation of noise between each unit is very low.Thus, the likelihood function of the map result is a product over all contributions, which is given by

Table 2 .
Performance of the proposed filter with different sizes of sub-area.