Optimal Particle Filter Weight for Bayesian Direct Position Estimation in a GNSS Receiver

Direct Position Estimation (DPE) is a rather new Global Navigation Satellite System (GNSS) technique to estimate the user position, velocity and time (PVT) directly from correlation values of the received GNSS signal with receiver internal replica signals. If combined with Bayesian nonlinear filters—like particle filters—the method allows for coping with multi-modal probability distributions and avoids the linearization step to convert correlation values into pseudoranges. The measurement update equation (particle weight update) is derived from a standard GNSS signal model, but we show that it cannot be used directly in a receiver implementation. The numerical evaluation of the formulas needs to be carried out in a logarithmic scale including various normalizations. Furthermore, the residual user range errors (coming from orbit, satellite clock, multipath or ionospheric errors) need to be included from the very beginning in the stochastic signal model. With these modifications, sensible probability functions can be derived from the GNSS multi-correlator values. The occurrence of multipath yields a natural widening of the probability density function. The approach is demonstrated with simulated and real-world Binary Phase Shift Keying signals with 1.023 MHz code rate (BPSK(1)) within the context of a real-time software based Bayesian DPE receiver.


Introduction
The Position, Velocity and Time (PVT) estimates of a Global Navigation Satellite System (GNSS) receiver are sometimes not sufficiently accurate, in particular in difficult environments such as urban areas, forests or indoors. In such environments, the observations do not follow a Gaussian distribution. This fact adds a bias to the PVT solution, which can be generally described as close to Gaussian distributed [1][2][3][4][5]. Thus, the PVT estimation in difficult environments can benefit from Sequential Monte Carlo (SMC) methods because they consider all available statistical information in the estimation process.
Particle filters [6,7] belong to the group of non-parametric Bayesian filters and are an implementation of an SMC method [8]. Thus, they allow for nonlinear and non-Gaussian distributed system and observation models. The performance comparison between particle filters and the Extended Kalman Filter (EKF) within a nonlinear and non-Gaussian environment shows that improved position estimates can be expected using particle filters with the drawback of computational intensive processing [9].
The application of particle filters in the field of PVT estimation can mostly be found in research. Often, particle filters are applied in tracking multiple objects (see e.g., [10]). For this purpose, particle  (1) refers to the analogue to digital conversion, the first stage of the receiver; (2) refers to the GNSS signal processing, which produces correlation values for each tracked channel as symbolically shown in (3). The correlation values are mapped and weighted to a particle cloud of a particle filter as shown in (4).
The first block (1) refers to the first stage of the GNSS signal processing. The Radio Frequency (RF) signal is gathered with a GNSS antenna. Afterwards, the RF signal is down-converted to an Intermediate Frequency (IF). Then, the analogue IF signal is converted into a digital signal by an Analogue to Digital Converter (ADC). Finally, the digital and to IF down-converted sample stream is transmitted via USB to a PC, where the whole GNSS receiver processing is performed.
The second step in the processing sequence is the signal correlation, indicated with (2) in Figure 1. At this stage, the receiver is configured to act in the same way as a vector tracking receiver. For each of the N GNSS signals, a complex-valued replica signal is generated using the best estimatef d and τ, at the Doppler and code phase. These estimates come from the latest PVT estimate (cf. step 4) after applying a LOS projection. Thus, there are Doppler and code phase estimates for each of the N channels. The signal replica is then correlated with the incoming signal to produce a new time series of correlation values. This is done for a wide range of code phase offsets, not only for Early, Prompt and Late, as it is typically done in GNSS receivers. Thereafter, a post-correlation Fast Fourier Transform (FFT) is applied to this complex valued correlation values in order to obtain synthetically generated correlation values in the Doppler direction. The resulting grid of correlation values in code phase and Doppler direction are called synthetically generated Multi-Correlation (MC) values (cf. step 3). These maps contain now correlation values for a range of Doppler bins and code phase offsets. The MC processing performs a navigation data bit wipeoff which also allows for long coherent integration times. The MC maps are centered at the old fed back Dopplerf d and code phaseτ, symbolically indicated with the black square. In the next step, each green PVT particle of a particle filter in the PVT domain (4) is mapped to its corresponding correlation value in (3). This is done for each channel N and using a sinc interpolation in both directions, code phase and Doppler. Based on the fact that the correlation values are coupled via the geometric relationship of its PVT particle, it is allowed to sum them up, the basic idea of DPE. The correlation values are used as input to determine the weight change of a particle, a crucial step in a Bayesian filtering approach, which is discussed in detail within this work. Using synthetically generated correlation values as input to the Bayes filter is a valid approach for DPE because it was proven in [26] that these synthetically generated correlation values represent mathematically the exact value of a real correlator. This approach mimics the direct correlation and it is computationally very efficient.
The last block (4) in Figure 1 refers to the navigation processor and thus to the implemented particle filter. There exist various types of particle filters-for example, Sampling Importance Resampling (SIR), Auxiliary SIR (ASIR) or Regularized Particle Filter (RPF). The details of these filters are not further discussed in this work, but all of them rely on a proper particle weight update. The particle weight update function determines the weight change of each PVT particle and takes the MC values as input. Afterwards, the final PVT solution (yellow circle in 4) is computed by constructing a probability density function using the particle weights and calculating the final PVT solution using a weighted mean over all particles. To illustrate the implementation, Figure 2 depicts the normalized logarithmic particle weights w i k from Equation (35) using the optimal particle weight update from Equation (34) (35), assuming a uniform distribution from the previous epoch, i.e., w i k−1 = 0 and using the optimal particle weight update from Equation (34) are shown. The particles are equidistantly distributed over a grid in the northeast plane. The lines through the plot correspond to the weighted correlation function in the position domain and thus refer to a GNSS signal. In a proper case (correct user velocity, clock error and drift), the lines overlap at a distinct point in the position domain, which is in this case the northeast plane. The resulting peak represents the probability of the 2D position. Note that the plotted weights are normalized and in the logarithmic scale, thus the peak has the maximum value of 0. The coherent integration time for this plot was set to T coh = 2 ms. The processed data refers to open sky. It was recorded at latitude LAT = 47.06446263 deg, longitude LON = 15.40777110 deg on the rooftop of Reininghausstraße 13a, Graz, Austria.
In BDPE, the implemented Bayesian filter plays an important role as it estimates the PVT solution from all available measurements. The following chapter discusses in detail the Bayesian framework and its role in the concept of DPE. . This plot depicts a case when clock error is improperly aligned. The GNSS signals do not overlap at a distinct position in the northeast plane. The plot shows the normalized logarithmic particle weights w i k from Equation (35), assuming a uniform distribution from the previous epoch, i.e., w i k−1 = 0 and using the optimal particle weight from Equation (34). The particles are equidistantly distributed over a grid. The lines through the plot correspond to the weighted correlation function in the position domain and thus refer to a distinct GNSS signal. The lines look very broad even after the weight update, which comes from the logarithmic scale given by log(w i k ). The coherent integration time for this plot was set to T coh = 2 ms. The processed data refers to open sky. It was recorded at LAT = 47.06446263 deg, LON = 15.40777110 deg on the rooftop of Reininghausstraße 13a, Graz, Austria.

Particle Filter
This chapter discusses briefly the basic principles of a particle filter in order to cover the complete framework of the investigated weight update step, which follows in the next chapters. Particle filters belong to the group of non-parametric Bayesian filters [13] and thus allow for dealing with nonlinear and non-Gaussian models [8]. In order to estimate the state of a dynamic system, two filter steps are required, Prediction and Update. The prediction step is based on the Chapman-Kolmogorov equation, which propagates the previous state x k−1 to the prior state x k by using the system model where z is a set of all available measurements, and f k (·) is a possibly nonlinear function with the corresponding process noise v k−1 , at epoch k. The process noise can be interpreted as the uncertainty of the system model, which is accounted for when propagating the system state. The incorporation of new information (measurements) is performed in the update step, which is described by the Bayes rule with the normalization constant The update incorporates the new measurement z k by using the measurement model, which relates the measurement z k to the states x k as where h k (·) is a possibly nonlinear function and where n k denotes the measurement noise, which describes the uncertainty of the measurement. Now, the basic idea of a particle filter is to recursively approximate the posterior probability function. If considering only the case that a filtered estimate is required at each time step, a simplified discrete approximation of the posterior probability density function (PDF) without dependency on the complete history can be written as where the posterior probability only depends on the previous state at epoch k − 1 and current measurement at epoch k. The approximation of the PDF in (6) is based on a set of support points {x i k , i = 1, ..., N s } with associated weights w i k , where the weights are normalized such that ∑ N s i=1 w i k = 1. The δ(·) in (6) denotes the Dirac delta measure. It is a crucial step of the particle filter to obtain the current epoch weight w i k , which depends proportionally on the previous epoch's weight w i k−1 which is propagated in time using by p(x i k |x i k−1 ) through (2), updated with new information by p(z k |x i k ) using the relation in (5) and normalized by the the importance density q(x i k |x i k−1 , z k ). The importance density q(·) defines the probability distribution of the support points, which is chosen in the presented implementation to be Gaussian distributed with an approximate mean of the truth in the very first step. The choice of q(·) is an important design step of a particle filter and one common approach is to choose the prior PDF such as which leads to the simple weight update equation The realization of a particle filter using (9) is known as a Bootstrap filter. The concept of a Bootstrap filter is also implemented within this work. To be able to update the weights optimally, a probabilistic measurement model for p(z k |x i k ) in (9) with dependency on the prior state must be available. In the proposed approach, the state vector x k at current epoch k is an 8-element vector containing the GNSS receiver position, velocity and time, where the position and velocity take three elements each and time two elements refer to the receiver clock error and receiver clock drift. The measurements z k at current epoch k refer to a vector containing the raw GNSS samples. In order to update the particle weight, BDPE directly maps the PVT at the particles to the raw signal samples, which contains a superposition of all GNSS signals. Thus, a probabilistic description of the samples dependent on multiple GNSS signals must be available. Such a probabilistic model exists in literature (see, e.g., (5.9) of [27]) for a single GNSS signal and is given as where the vector s ∈ C L contains all signal samples s µ , µ ∈ {1, ..., L}, which depend on the signal model parameters, the real valued signal amplitude A ∈ R + , the time delay τ ∈ R, the Doppler frequency ω ∈ R , the uniformly distributed carrier phase φ ∈ [0, 2π[ and i the imaginary number. The term c(t µ − τ) ∈ [−1, 1] refers to the PRN code and the term exp{iωt µ − iφ} to the signal carrier. The parameter t µ ∈ R defines the time at sample index µ and depends on the sampling frequency f s ∈ R such as t µ = µ/ f s . The probability distribution p(z k |x i k ) in the weight update Equation (9) can be related to the sample distribution in (10), where z k equals the measured sample vector such as s = z k , where t µ=1 corresponds to the time of the first sample in epoch k and t µ=L to the time of the last sample in epoch k. The parameters A, τ, ω, φ for a single replica can directly be related to the PVT state z k , respectively, to the state of each PVT particle. The state z k is defined as where the three-dimensional position r rx and three-dimensional velocity v rx relate to the Earth Centered Earth Fixed (ECEF) reference frame. The time vector t rx contains the GNSS receiver clock error t rxClkErr and clock drift t rxClkDri f t . The replica signal code phase τ n and Doppler ω n for a single signal of satellite n are related to z k using τ n = geom. distance ||r sat,n − r rx || c + t rxClkErr + T n + I n , where r sat,n and v sat,n correspond to the position and velocity of satellite n, c to the speed of light, T n to the tropospheric delay, I n to the satellite and frequency dependent ionospheric delay and f c to the carrier frequency. These relationships assume that the receiver position is known at epoch k and thus the ephemeris of the satellite n for the generated replica signal needs to be available. In a simplified simulated setup without ionosphere, troposphere and other errors, the terms I n and T n may reduce to zero and can be neglected. The simulations which are performed within this work do not consider any geometric relationship, troposphere and ionosphere for analysis, but the real-world examples fully consider the geometric relationship as well as the ionospheric and tropospheric terms. Basically the amplitude A and carrier phase φ can also be estimated from the particle filter, but this further increases the dimensionality of the estimation problem form currently 8 to 10 dimensions. To avoid a further increase in dimensionality and thus computational complexity and be able to take the synthetically generated multi-correlation values which only depend on code-phase τ and Doppler ω, the proposed approach is to make the probability function independent of the signal amplitude A and carrier phase φ, which is discussed in the next chapter.

Derivation of the Optimal Particle Weight for Multiple GNSS Signals
Basically, the goal of the derivation is to retrive a formulation for an optimal particle weight in dependence of an abitrary number of used GNSS signals dependent on the code phase τ and Doppler ω only. Thus, the following equations show the extension of the probabilistic model from single to multiple GNSS signals in Equation (17) with dependency on A, τ, ω, φ, in Equation (24) with dependency on A, τ, ω and in Equation (27) with dependency on τ, ω.
Under the assumption that the carrier phase is uniformly distributed over 0 and 2π, it can be integrated out as shown in [27] to with the definition of the correlator value P n (τ, ω) as where the correlation value P n (τ, ω) is taken from the multi-correlation maps as shown in Figure 1 and I 0 referes to the Bessel function of the first kind and order zero. A similar Equation for joint detection of weak GNSS signals can be found in [18]. In order to extend the model to multiple GNSS signals, let's do this in the first step for two signals and find later a mathematical rule to extend it to multiple signals. The basic model for two signals is defined as and with the definition of |z| = √ zz, |z| 2 = zz and the redefinition of the replica signals r 1 and r 2 , m can be expanded to With |r n | 2 = A 2 n |c(t µ − τ n )| 2 and neglecting the cross-correlation terms leads to Insert Equation (21) in Equation (17), consider that the sum of the squared PRN code leads to ∑ L µ=1 |c(t µ − τ 2 )| 2 = L and realigning the terms to integrate out φ 1 and φ 2 leads to It can be shown in Equation (24) that the extension to multiple GNSS follows the same structure as shown in Equation (14) for a single GNSS signal, if the orthogonal cross correlation terms are neglected as shown in Equation (20).

Solving an integral of type
c }I 0 (xp)dx was performed with the help of Wolfram Mathematica (Version 11.3, Wolfram Research Inc., Champaign, IL, USA) [29], which helps to integrate out the amplitude as With the assumption of using Gaussian distributed complex valued samples s µ whose real and imaginary parts are each of variance one and zero mean, the sum over the absolute squared samples approximates to L ≈ 1 2 ∑ L µ=1 |s µ | 2 . This approximation can be made if the GNSS signal amplitude is significantly smaller than the noise. This is generally the case for GNSS signals when considering receivers operated on the ground [1]. With this approximation, Equation (26) can be rewritten for N GNSS signals as Equation (27) defines the probability function p(z k |x i k ) for an arbitrary number N GNSS signals in order to obtain w i k in Equation (9). Note that the code-phase and Dopplers changed now to vectors τ ∈ R N and ω ∈ R N . However, with this equation, the weight can not be evaluated directly due to finite precision effects. In particular, the terms exp |P n (τ n ,ω n )| 2 4 and I 0 |P n (τ n ,ω n )| 2 4 use a typically 'large' correlation value P, which causes a numerical problem. The common maximum value that can be stored on 64-bit platforms is a double precision floating point value (not considering dedicated floating point libraries). The maximum value is limited to 1.7E ± 308 (15 digits) [30]. Considering the term exp |P n (τ n ,ω n )| 2 4 in Equation (27), the numerical limits are reached for a correlation value P n (τ n , ω n ) > 53.28 approximately. A solution to overcome the effect of limited digital precision is proposed in the next chapter.

Logarithmic Weight Update
When performing the particle weight update, the numerical values can be too large to be evaluated by a computer. The numerical issues at the weight update when evaluating Equation (27) can be solved by shifting the equations to the logarithmic scale and performing an additional normalization step. As a first step, let be an approximation of the Bessel function of first kind where only the first term is used to approximate the function with order zero as With Equation (29), Equation (27) can be rewritten to Based on the reason that a later introduced normalization step is performed in Equation (35), the constant term H const can be neglected. With that, the weight update Equation in (9) can be expressed as Let us define w i k = log(w i k ) and w i k−1 = log(w i k−1 ); then, the update can be written in logarithmic scale as The weight update from Equation (9) can now be done in five steps from Equation (33) to Equation (37) using This series of equations allows for performing the weight update within numerical boundaries (e.g., double precision) and can generally be used to implement the optimal particle weight update in any Bayesian filter working with correlator values.

Discussion of the Ranging Accuracy
A further investigation determined that the probability distribution after the weight update from Equation (36) shows a very small variance in range, especially for strong GNSS signals and long integration times (e.g., larger than ten milliseconds). This is briefly illustrated in Figure 4. It can be seen that the amplitude of the correlation value has a significant impact on the variance of the probability function. For a typical open sky case at 45 to 50 dB-Hz and integration times of ≥10 ms, a weighted standard deviation of ≤0.32 m would be observed. This is caused by the large amplitude of the correlation function and small variations in the range cause already a significant variation in P(τ).
This behavior is expected and consistent with standard tracking theory and for example represented by the Cramér-Rao lower bound (CRLB) of the code pseudorange noise. However, the ranging accuracy described by the CRLB is never reached in practice, as multipath, residual orbit or satellite clock errors or residual ionospheric delays by far exceed the CRLB. Least-squares positioning algorithms or a Kalman filter (both working with code pseudoranges) can cope well with this increased ranging errors as the overall ranging accuracy is assumed to be unknown and is estimated from the code pseudorange residuals. Those algorithms only consider the relative accuracy variations between the different pseudoranges to define the relative weights, but do not account for a common accuracy scale factor.
However, for DPE, it is expected that this narrow probability distribution of the weights cause problems. First, grid based filters with fixed resolution and extension may not be able to cover the narrow probability distribution appropriately. Second, a residual code delay bias (multipath, orbit, clock or iono) may shift the probability density of one satellite with respect to another. Instead of overlapping the true PVT, the shifted functions may assume very small values (or even be zero due to the finite precision in the computer) and thus eliminate particles at the true PVT. The particle cloud of the particle filter would be severely concentrated around one maximum and it is expected that the filter may also have difficulties to track the PVT solution accurately in case the amplitudes of the different GNSS signals vary in time relative to each other.
Thus, a trade-off and an adaptive approach for the grid resolution, grid extension and computational power seem to be necessary to guarantee a good handling of the narrow probability functions. One could think of increasing the process noise to cope with the delay biases, but as the process noise is determined by the dynamic model, this is not a viable approach.
The next chapter proposes a method to overcome these problems by introducing an unknown delay bias as nuisance parameter.

Delay Bias as Gaussian Nuisance Parameter
This section describes an approach to handle user range uncertainties (due to e.g., multipath, orbit, clock and iono) within a DPE framework, which are much larger than the code noise. The idea is to introduce an unknown normally distributed bias in the signal model which reflects uncertainties in the user range. This artificially degrades the accuracy of the raw signal but better reflects the true measurement situation and avoids situations of non-overlapping probability functions.
One may argue that this delay bias degrades the overall positioning accuracy, but this argument is only true for an ideal situation of vanishing range uncertainties. For real situations, Equation (27) does not reflect the range uncertainties and thus is therefore only an approximate model of the signal samples. Overall, we think that the parameter can be justified and its introduction produces more meaningful probability density function compared to the ideal case.
Basically, this bias can be applied at code phase τ and Doppler f d , referred to as code delay bias ∆τ and Doppler bias ∆ f d . This work discusses and analyses the introduction of a code delay bias ∆τ to the code delay τ → τ + ∆τ, but the procedure for the Doppler bias is basically the same. The code delay bias ∆τ should express orbit errors, ionospheric errors, clock errors and other modeling errors. It is assumed that these errors follow a normal distribution such that Considering for now only one signal and applying the delay bias from Equation (39) to Equation (26) leads to As this integral cannot be further simplified analytically, we target an approximate evaluation at selected grid points ∆τ k . Neglecting a normalization factor Equation (40) can be expressed as The grid points ∆τ k should be chosen in a way to cover a significant part of N(0, σ 2 ∆τ ) from Equation (38) with a reasonable resolution. In this work, we choose a range of ±3σ ∆τ with a resolution of 0.01 m. Under the assumption that the introduced code delay bias follows the same distribution for both signals, Equation (41) can be rewritten to Applying the same simplification steps as was done to obtain Equations (27) and (30) Equation (43) is also numerically difficult to evaluate due to the large exponential. The sum over the grid points does not allow directly shifting the Equation to the logarithmic domain in order to get rid of the exponential, thus some reformulation is necessary. Let us define Then, the inner term of the sum can be rewritten for K = 2 to 2 ∑ k=1 c k = a 1 e b 1 + a 2 e b 2 = a 1 e b 1 which can be written in a general form as When selecting b 1 to be the maximum, the exponential term in the sum of Equation (47) remains < 1. Applying the logarithm now leads to which can be well evaluated numerically. Using this in Equation (43) and applying the logarithm leads to the new logarithmic weight update equation Equation (50) can now be used instead of Equation (35) within the weight update sequence Equation (33) to Equation (37) to include the Gaussian delay bias. An additional sum over K grid points and a maximum search over b k needs to be evaluated for each signal N, which increases the computational complexity. It should, however, be noted that the computation needs to be performed only for each set of multicorrelator values for each satellite signal and then an interpolation can be performed to obtain the weight update for the individual particles. This avoids computation of Equation (50) for each particle.

Simulations and Real-World Results
The previous chapters presented the optimal weight update function in Equation (32) for DPE when using a Bayes filter and dicussed in Figure 4 the problem with the resulting narrow probability function. In Equation (50), a possible solution is proposed to overcome the presented issues. In order to understand the behaviour of the proposed weight update function, an analysis is performed for Equation (50) within this chapter for a single signal N = 1 and assuming a uniform distribution from the previous epoch, i.e., w i k−1 = 0. All simulations have been performed completely with MATLAB, except the real-world scenarios, which take as input the correlation values from the GNSS software receiver. The MATLAB simulation performs a correlation with a simulated GNSS signal. In order to see clearly the impact of the different parameters, the simulated signal was generated without noise. All following plots show the correlation values |P| from Equation (15) and the resulting weight w i k from Equation (36) when using Equation (50) as weight update function. For visualisation, Equation (36) is plotted instead of Equation (37). This is because of the normalized amplitude to 1, which makes a comparison more easy to visualize. It is noted that the analysis is based on the herein presented equations. The derived equations do not consider any present multipath signal, thus there is no estimation of multipath parameters nor a handling or mitigation of multipath.
The next subsections analyse the impact on the resulting probability function for

Impact of Different Signal Strengths
Basically, a higher GNSS signal amplitude increases the correlation value from Equation (15) and thus influences Equation (50). From this, it is expected that changes in amplitude do not shift the mean of the resulting probability function but influence the variance. This is obvious because a stronger received GNSS signal must result in a more accurate estimate, which is shown in Figure 5.

Impact of Different Code Delay Bias Variances
Basically, the expected measurement accuracy σ w strongly depends on the amplitude of the correlation value |P| and thus on the coherent integration time T coh and the GNSS signal strength. However, with the introduction of a theoretical lower limit of the measurement accuracy in Equation (40) defined by σ ∆τ , the estimated accuracy σ w must approach the lower limit σ ∆τ with an increasing correlation value |P|. This behaviour is shown in three consecutive plots in Figures 6-8 with an increasing coherent integration time T coh = [1, 10, 100] ms. From this series, it is clearly visible that σ w also approaches a very small σ ∆τ = 0.1 m in the case of a long coherent integration time T coh = 100 ms.
The magnitude of σ ∆τ should cover residual user range errors (coming from orbit, satellite clock, multipath or ionospheric errors), which can be in the range of several meters. . Impact of different code delay bias standard deviations σ ∆τ on the weights w i k for a coherent integration time T coh = 1 ms. For the given conditions, only the weighted standard deviation σ w for σ ∆τ = 6 seems to approach the theoretical lower limit of σ ∆τ . In particular, σ ∆τ = 0.1 can not be reached due to the influence of the low correlation time T coh at given C/N0.   Figure 8. Impact of different code delay bias standard devications σ ∆τ on the weights w i k for a coherent integration time T coh = 100 ms. It can be seen that σ w now approaches the theoretical limits for all σ ∆τ . Note also the significantly changed amplitude on |P|.

Impact of Constructive and Destructive Multipath
The influence of constructive and destructive multipath is shown for a realistic case with a multipath amplitude of α MP = 0.5 with respect to the LOS signal and a multipath offset of ∆τ MP = 50 m in Figure 9. As similar to the conventional Delay-Lock-Loop (DLL) in a GNSS receiver and thereof resulting pseudorange measurement, the weighted mean µ of the probability function also becomes biased. It can be seen that the resulting bias in µ is in the opposite direction for the constructive and destructive case. Furthermore, the destructive multipath reduces the amplitude of the correlation value |P|, which increases the weighted standard deviation σ w , while the constructive multipath acts in the opposite direction and may also reduce the variance in distinct cases as visible in Figure 10 for the black line. Furthermore, it can be seen that constructive and destructive multipath changes the resulting probability function in a different way, even if the multipath parameters ∆τ MP and α MP are exactly the same.

Impact of Short, Medium and Far Multipath
The impact on short, medium and far multipath on the probability function is shown on a strong multipath case α MP = 0.9. Basically, it is expected that such a strong multipath seldomly occurs in typical real-world scenarios, but it was chosen to clearly visualize the impact of an increasing multipath offset. It can be seen in Figure 10 that an increase in ∆τ MP directly increases σ w .
An interesting effect is the short constructive multipath case with an ∆τ MP = 3 m. In this case, the weighted standard deviation σ w is underestimated and shows a smaller variance than the LOS signal. This can be dangerous because this leads to a biased and at the same time more accurate measurement, compared to the truth. Considering the present multipath and incorporating it into the models for the weighting may prevent such an underestimation.   Figure 10. Impact of different multipath offsets ∆τ MP on the probability function. For a better visualization, a strong relative multipath to the LOS signal with α MP = 0.9 is chosen. It can be seen that higher offsets increase the weighted mean µ but do not necessarily increase the weighted standard deviation σ w . The dotted black line refers to the LOS signal.

Impact of Different Multipath Amplitudes
It is expected that an increase of the multipath amplitude α MP further shifts the mean and increases the variance. This case is outlined in Figure 11. Naturally, in the case that the multipath signal is as strong as the LOS signal, the resulting probability function automatically covers the whole uncertainty range with a mean exactly between the two signals. This case is shown with the dash-dotted blue line in Figure 11.  Figure 11. An increase of the multipath amplitude α MP shifts the weighted mean µ and increases the weighted standard deviation σ w . The probability function naturally covers the uncertainty also in the case of α MP = 1, when the multipath signal is as strong as the LOS signal. The dotted black line refers to the LOS signal.

Real-World Open Sky and Urban Scenario
The real-world tests evaluate the probability function with realistic data. Therefore, two scenarios have been selected: (1) an open sky case on a rooftop antenna in order to see the probability function under ideal conditions and (2) an urban case, where multipath is present and a significant shift and deformation on the probability function should be present.
The data was recorded in Graz, Austria and post-processed with the SX3 software receiver. The receiver was configured to dump the multi-correlator maps to files, which contain the correlation values |P|, which are related to the latest PVT. The multi-correlator maps contain correlation values for a range of code phases τ and Dopplers f d , as shown in (3) of Figure 1. For all real-world tests, the red crosses in the upper plots correspond to the correlation values along the Doppler bin, which contains the correlation maxima. In order to obtain a high-resolution probability function, the correlation values (red crosses) have been interpolated with sinc function (black line).
It is of major importance that the input correlation values |P| are of correct amplitude because it basically drives the shape of the probability function and variance estimate. In order to obtain the same amplitude for |P| as in Equation (15), all amplitude scaling elements have been verified in the receiver. The GNSS Signal Sampless µ in Equation (15) have been recorded with a quantization of 2 bits in the Analog to Digital Conversion (ADC) stage. The two bits refer to the value range [−3 −1 1 3] and thus the RMS of the recorded samples differs from that in Equation (15) because the Automatic Gain Control (AGC) within the receiver front-end steers the amplitude of the GNSS signal in a way to optimally use the available quantization range. The RMS value of the sampless µ can be measured by the receiver. In this experiment, the RMS for the GPS L1 band was measured with β s µ ,RMS = 1.71, as also shown in Figure 12. Due to computational efficiency, the replica signal used for the correlation was generated with an amplitude of 8, thus it is simply β rep = 8. Based on this, the correlation values can be obtained after applying the scaling factors as where |P MC,map (τ, ω)| refers to the correlation value from the multi-correlator map of the software receiver. The urban scenario was recorded during a test drive through Graz. The data has been post-processed and analyzed with the software receiver. An urban environment with significant deformation on the correlation function was selected and analyzed with MATLAB. The analysis focuses on the GPS L1 C/A satellites PRN 27, PRN 21 and PRN 18, as respectively shown in Figure 14, Figures 15 and 16. A overview of the analyzed scenario is given in Figure 17, which shows the three analyzed positions (red dots) and the satellite constellation. Satellite PRN 27 in Figure 14 and satellite PRN 18 in Figure 15 show a significant variation in amplitude of the correlation values |P|, basically due to shadowing and multipath effects. Satellite PRN 21 in Figure 15 is less affected by the environment, because the signal amplitude remains nearly constant and no other significant effects are visible. Nearly all of the probability functions face a significant bias in µ and increased weighted standard deviation σ w compared to a standard open sky signal. Especially the weighted standard deviation of PRN 18 at third column in Figure 16 shows a significant increase due to the lowered correlation amplitude. Furthermore, the correlation function as well as the probability function of satellite PRN 27 in Figure 14, position at time W/S 1901/317840.4 (middle row) seem to be significantly affected by multipath, which is assumed due to the flattened peak in the correlation function and the resulting deformation of the probability function.   Figure 14. This satellite signal is significantly affected by the environment. From the azimuth of satellite PRN 18 and the location of the buildings as shown in Figure 17, it can be assumed that the GNSS signal is blocked at the first and last positions, which fits to the amplitude of the correlation values. Interestingly, in the case of the third position, the significantly small correlation value leads to an increase in variance estimate.
For clarification and as noted briefly above, all presented real-world scenarios show the resulting probability function with respect to the latest PVT estimate, not to an high accurate absolute reference.
It should be considered that the latest PVT estimate might be already biased and thus the herein presented offset µ might not show the true bias. Dedicated analysis is planned for future studies using in one case a GNSS simulator and in another case a surveyed static reference position together with precise orbits. Nevertheless, the presented figures give an impression of the effects occurring (deformations and shifts) on the probability function. This is clearly visible when comparing the open sky case to the urban scenario, even if the biases can not be quantified with absolute values.

Conclusions
This paper discusses Bayesian Direct Position Estimation and describes a real-time capable implementation in a software based receiver. The Bayesian estimation framework is described in detail and a link to direct position estimation is established. The work shows the derivation of the optimal particle weight for BDPE and delivers a solution for performing the weight update in the logarithmic scale to overcome limited precision of digital computation. It further points out problems within Bayesian filters when performing the optimal particle weight update directly. Therefore, the paper proposes a method that introduces an additional nuisance parameter on the code measurement to cover the residual user range errors. The impact on several scenarios for a single GPS C/A BPSK(1) signal is given. These scenarios include a brief discussion on the impact of different signal strengths, different code delay bias variances, constructive and destructive multipath and different multipath amplitudes. Additionally, two real-world examples show the results of the optimal particle weighting applied on ideal open sky data and on multipath-prone urban environment data.
The work extends the original particle filter from [23], which followed a heuristic approach for the particle filter weight update. The original filter used a suboptimal chosen exponential weight update function, which causes a dense concentration of the particle cloud. In order to avoid such a dense concentration, the process noise was increased to cover also the residual user range errors. This is a working and viable approach to cover these errors, but it has caused rapid resampling steps of the particle filter.
From the proposed method, it is expected that a Bayesian filter can perform an optimal particle weight update with numerical stability. Furthermore, it is expected that the particle filter positioning performance and robustness will improve. The positioning benefits from the optimal particle weight update and it is expected that the estimated position becomes more accurate and the estimated accuracy of the position becomes more reliable. The introduced nuisance parameter allows for overlapping GNSS signals in the PVT domain in the case of unmodeled residual user range errors, which is expected to make, on one hand, the filter more robust and, on the other hand, reduce the number of resampling steps.
Author Contributions: J.D. implemented the particle filter, derived the formulas, did the analysis and wrote the paper. K.F. supported the mathematical derivation of the equations. T.P. provided the concept for the introduced nuisance parameter and initiated the work.
Funding: This research received no external funding.

Acknowledgments:
The authors would like to thank IFEN GmbH, Poing, Germany, for providing the GNSS Software Receiver. Moreover, we want to acknowledge Professor Manfred Wieser from the Graz University of Technology, Institute of Geodesy, Austria, for helpful discussions on BDPE.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: