1. INTRODUCTION
Global Positioning System (GPS) technology plays an important role in positioning and navigation applications due to the increasing demand for a positioning system that can provide continuous, accurate and reliable position. However, this may not always be possible due to different environmental effects on GPS signals such as attenuation, fading and multipath (Kaplan and Hegarty, Reference Kaplan and Hegarty2005). Therefore, efforts need to be made to make the civilian part of the system more accurate, reliable and available for urban canyon and indoor navigation environments in particular.
Despite the major advances in signal processing methods and technologies used nowadays inside GPS receivers, multipath remains a major source of ranging error (Soubielle et al., Reference Soubielle, Fijalkow, Duvaut and Bibaut2002). The delay of multipath signal depends on the extra time of the travelled path by the reflected signal, which is longer than the direct signal's path. In deep urban canyons, high-rise buildings and surrounding objects near the receiver lead to the reception of indirect (reflected) signals in addition to the direct Line of Sight (LOS) signal. This effect can significantly distort the shape of the correlation function used for Time-Of-Arrival (TOA) estimation inside the GPS receiver (Broumandan and Lin, Reference Broumandan and Lin2008). Due to inaccurate LOS delay estimation in multipath environments, receivers cannot provide reliable Position, Velocity and Time (PVT) solutions.
One specific circumstance in multipath propagation is the situation when the LOS component overlaps with one or several short-delay multipath components making the TOA estimation process more complicated (Sleewaegen and Boon, Reference Sleewaegen and Boon2001). This short-delay multipath scenario commonly happens in urban environments and is the main challenge addressed in this paper.
The most familiar code tracking algorithms used inside commercial GPS receivers are the classical delay tracking techniques. These techniques include the feedback delay estimator loop known as the Delay Lock Loop (DLL) and the use of special correlators and discriminator functions. The DLL estimates the code delay by maximising the correlation between the incoming signal and the local code replica (Parkinson and Spilker, Reference Parkinson and Spilker1996). The structure of the DLL consists of one or more correlators, a non-linear discriminator for evaluating the code delay error function, a loop filter for improving code delay estimation and a code numerically controlled oscillator (Kaplan and Hegarty, Reference Kaplan and Hegarty2005).
In the standard DLL discriminator function, two correlators called Early (E) and Late (L) correlators are spaced at a code chip width from each other and they are used to design the discriminator function. The code delay is estimated by balancing the amplitude/power on the E and L correlators. This method is called Early-Minus-Late (EML) technique (Kaplan and Hegarty, Reference Kaplan and Hegarty2005). The classical EML fails to mitigate the multipath effects. Thus, several enhanced EML-based techniques are introduced in the literature (Bhuiyan and Lohan, Reference Bhuiyan and Lohan2010) to reduce the multipath effects, especially the short delay multipath scenarios.
The Multipath Estimation Delay Lock Loop (MEDLL) is another class of code tracking algorithm (Van Nee and Siereveld, Reference Van Nee and Siereveld1993). It is a parametric mitigation technique that attempts to minimise the effect of multipath by estimating the parameters that characterize the received signal. NovAtel Inc. implements it for GPS receivers. MEDLL aims at jointly estimating the amplitudes, delays and phases of both the LOS component and multipath components. MEDLL uses many correlators to accurately determine the shape of the distorted correlation function (Townsend et al., Reference Townsend, Fenton, Van Dierendonck and Van Nee1995). The main idea is that the multipath components are estimated recursively and the effect of each path is removed before the estimation of the next. Several modified MEDLL techniques have been proposed in the literature to reduce the computational cost (Sokhandan, Reference Sokhandan2013). The MEDLL technique outperforms the classical multipath mitigation techniques but it does not completely eliminate all multipath errors (Bhuiyan et al., Reference Bhuiyan, Lohan and Renfors2008). Multipath signals with short relative delays, specifically, are difficult to reject.
The ultimate goal of this research is to develop an innovative GPS tracking technique to mitigate multipath effects in challenging environments such as in urban canyons where GPS receivers suffer from closely spaced multipath signals. The main objective of the research is to develop a new approach for multipath mitigation based on Fast Orthogonal Search (FOS), offering a higher resolution code tracking performance than the conventional techniques in harsh environments.
FOS is utilised in this research to model the distorted multipath correlation function and obtain better TOA estimation for both LOS and multipath signals for GPS receivers. Gram–Schmidt orthogonalisation is used to create a series of orthogonal functions from a given set of non-orthogonal candidate functions. The candidate functions are designed in this research to represent the GPS L1 reference correlation functions with several multipath delays, phases and amplitudes. The proposed method tries to find a perfect match between the distorted correlation function and a wider selection of FOS pre-designated candidate functions.
2. BACKGROUND ON MULTIPATH MITIGATION TECHNIQUES
2.1. Signal Model
The GPS signal at the output of the GPS receiver antenna can be written as (Borio, Reference Borio2008):
$$r\left( t \right) = \mathop \sum \limits_{i = 1}^N {A_i}{d_i}\left( {t - \tau _i^0} \right)\; {c_i}\left( {t - \tau _i^0} \right)\; exp\left\{ {\,j2\pi \left( {{\,f_L} + f_{D,i}^0} \right)t + \phi _i^0 (t)} \right\} + MP\left( t \right) + \eta (t)$$
where r(t) is the received signal that is the sum of N useful signals that are broadcast by N different GPS satellites, A i is the received carrier power including the effect of atmospheric attenuation and receiver antenna gain pattern, d i(t) represents the data navigation massage, 
$\tau _i^0 $ is the LOS path delay introduced by the channel, c i(t) represents the product of Pseudorandom binary sequence (PRN) code, secondary code and sub-carrier and f L is the carrier frequency which depends on the GPS band under consideration. For GPS L1 band, f L = 1575·42 MHz and for GPS L2 band, f L = 1227·60 MHz and for GPS L5 band, f L = 1176·45 MHz. 
$f_{D,i}^0 $ is the Doppler frequency affecting the ith useful signal, 
$\phi _i^0 (t)$ is a time varying phase offset, MP(t) represents the sum total of multipath reflections and η(t) is an Additive White Gaussian Noise (AWGN) with single sided noise Power Spectral Density (PSD) of 
$${{N_0}/ 2}$$ W/Hz.
After the GPS signals reach the receiving antenna, the received signals are filtered, down converted to baseband and digitised by the Radio-Frequency (RF) front-end (Borre et al., Reference Borre, Akos, Bertelsen, Rinder and Jensen2007). Neglecting the quantisation effect, the GPS signal at the output of the RF front-end can be rewritten as:
$$\eqalign{r\left[ n \right] = \mathop \sum \limits_{i = 1}^N {A_i}{d_i}\left[ {n - \tau _i^0} \right]{c_i}\left[ {n - \tau _i^0} \right]exp\left\{ {\,j2\pi \left( {{\,f_L} + f_{D,i}^0} \right)n{T_s} + \phi _i^0 \left[ n \right]} \right\} + MP\left[ n \right] + \eta \left[ n \right]}$$
where r[n] = r[nT s] is a discrete-time sequence r[n], obtained by sampling a continuous-time signal r(t) with a sampling frequency 
${{f_s} =1/T_s}$
The delay tracking is based on evaluating and processing the correlation 
${\rm {\cal R}}(.)$ between the incoming signal and each PRN code across all possible combinations of local code offset and Doppler shift and can be defined as (Borio, Reference Borio2008):
$${{\rm {\cal R}}_c}\left( {\bar \tau, {\rm \;} \overline {{\,f_D}}, n} \right) = {\rm \;} \displaystyle{1 \over {{N_c}}}{\rm \;} \mathop \sum \limits_{t = 0}^{{N_c} - 1} r\left( {n - t} \right){\rm \;} c\left( {n - \bar \tau - t} \right){\rm \;} exp\left\{ {\,j2\pi ({\,f_{IF}} + \overline {{\,f_D}} )(n - t){T_s}} \right\}$$
where the parameter N c is the number of samples used for computing a single correlation output and T c = N cT s is the coherent integration time in milliseconds for the GPS signal. The correlation is taken by using the last N c samples until the current time index n. 
$c\left( {n - \bar \tau - t} \right)$ is the local replica reproducing the C/A code, 
$\bar \tau $ and 
$\overline {{f_D}} $ are the code delay offset and the Doppler shift, respectively, tested by the receiver. In order to handle the noise, non-coherent integration can be used (Yang et al., Reference Yang, Miller, Blasch and Nguyen2001). The non-coherent average correlation function 
${{\rm {\cal R}}_{NC}}\left( {\hat \tau, \; {{\hat f}_D}} \right)$ can be written as:
$${{\rm {\cal R}}_{NC}}\left( {\hat \tau, \; {{\hat f}_D},n} \right) = \; \displaystyle{1 \over M}\; \mathop \sum \limits_{m = 0}^{M - 1} {\left \vert {{{\rm {\cal R}}_c}\left( {\bar \tau, \; \overline {{\,f_D}}, n - m} \right)} \right \vert ^2}$$where M is the total number of correlation samples that are non-coherently integrated. It should be noted that for M = 1 only coherent integration is used.
In a classic GPS receiver, three replicas for each PRN code are generated named prompt, early, and late. The resulting correlation outputs of prompt (p) replica of the PRN code in the presence of number of L multipath components are written as (Bhuiyan et al., Reference Bhuiyan, Lohan and Renfors2008):
$${I_{\,p,n}} = {A_n}{d_n}\left( {\tau _n^0} \right){\rm {\cal R}}\left( {\tau _n^0} \right){\it cos} \phi _n^0 + {A_n}\mathop \sum \limits_{l = 1}^L {\alpha _l}{d_n}\left( {\tau _n^0 - \Delta {\tau _l}} \right){\rm {\cal R}}\left( {\tau _n^0 - \Delta {\tau _l}} \right)cos(\phi _n^0 + \Delta {\phi _l})$$
$$\eqalign{Q_{p,n} & = \; {A_n}\; {d_n}\left( {\tau _n^0} \right){\rm \; {\cal R}}\left( {\tau _n^0} \right)\; Sin\; \phi _n^0 + \; {A_n}\; \mathop \sum \limits_{l = 1}^L {\alpha _l}\; \; {d_n}\left( {\tau _n^0 - \Delta {\tau _l}} \right){\rm \; {\cal R}}\left( {\tau _n^0 - \Delta {\tau _l}} \right)\; sin({\phi}_n^0 \cr & \quad + \Delta {\phi}_l)}$$where α l is the relative amplitude of l-th multipath, Δτ l is the relative delay and Δϕ l is the relative phase with respect to LOS signal, p refers to prompt correlator, n is the current time index and L is the multipath components.
2.2. Early-Late-Slope (ELS) technique
The basic idea of the ELS technique is to determine the slopes of the early and the late sides of the correlation function's central peak (Brodin and Daly, Reference Brodin and Daly1997). Based on the slope values, the code-phase correction term is calculated, which goes back to the DLL algorithm to shift the pairs of correlators. In the existence of multipath, the slopes at the early and late sides of the autocorrelation function of the received signal are no longer of the same magnitude. Thus the code-phase correction term is computed in order to compensate the slope difference of the rising and falling edges of the distorted autocorrelation function and it can be expressed as (Bhuiyan et al., Reference Bhuiyan, Lohan and Renfors2008):
$$\Delta \tau = \displaystyle{{{y_1} - {y_2} + d/2({a_1} + {a_2})} \over {{a_1} - {a_2}}}$$
where Δτ is the code tracking error (Δτ will equal zero when the two correlators are located at equal distance on each side of the correlation peak), 
${a_1} = \displaystyle{{2({E_1} - {E_2})} \over d}$ and 
${a_2} = \displaystyle{{2({L_1} - {L_2})} \over d}$ are the autocorrelation function slopes at the early and late side, respectively, y 1 and y 2 are the values of autocorrelation function at E 1 and L 1, and d is the spacing between E 1 and L 1. The spacing between E 2 and L 2 is equal to double the spacing between E 1 and L 1.
2.3. High Resolution Correlator (HRC)
The High Resolution Correlator (HRC) utilises the outputs of five correlators (E 1, E2, P, L1, L2) to form a linear combination of correlator outputs that give a net correlation function that is much narrower than the typical C/A code autocorrelation function (Irsigler et al., Reference Irsigler, Hein and Eissfeller2004). The E 2 and L 2 correlators are spaced twice the space of E 1 and L 1 correlators by 2d chips and the P (prompt) correlator central. Initially, a created correlator, named HRC Prompt, is modelled as follows (Sokhandan, Reference Sokhandan2013):
$${P_{HRC}} = 2P - ({E_1} + {L_1})$$Early and late HRC correlators are calculated to form an EML discriminator function for code tracking as follows:
$${E_{HRC}}\left( \tau \right) = {P_{HRC}}\left( {\tau + d} \right) = 2{E_1} - \left( {{E_2} + P} \right)$$
$${L_{HRC}}\left( \tau \right) = {P_{HRC}}\left( {\tau - d} \right) = 2{L_1} - \left( {{L_2} + P} \right)$$Finally, the corresponding EML discriminator function can be expressed as
$${D_{HRC}}\left( \tau \right) = \displaystyle{{\left( {{E_{HRC}} - {L_{HRC}}} \right)} \over 2}$$
$${D_{HRC}}\left( \tau \right) = \; \left( {{E_1} - {L_1}} \right) - \left( {{E_2} - {L_2}} \right)/2$$The HRC technique helps only to reduce code multipath. A modified version of HRC was proposed by (Garin and Rousseau, Reference Garin and Rousseau1997) to solve the issue of the reduction in signal power, which degrades the tracking accuracy at low C/N0 values (McGraw and Braasch, Reference McGraw and Braasch1999). This modified implementation uses the following discriminator function:
$$\Delta {\tau}_{HRC} = \displaystyle{{I{D_{HRC}} \times IP + \; Q{P_{HRC}} \times QP} \over {SP}}$$where ID HRC is the real part of D HRC, IP and QP is the real and imaginary parts of the prompt correlator, QP HRC denotes the imaginary part of P HRC and SP = IP 2 + QP 2. This implementation of HRC reduces carrier phase multipath as well as providing an enhanced reduction in code multipath.
2.4. Multipath Estimation Delay Lock Loop (MEDLL)
MEDLL is an advanced multipath mitigation technique, which was proposed by Townsend et al. (Reference Townsend, Fenton, Van Dierendonck and Van Nee1995) and Van Nee and Siereveld (Reference Van Nee and Siereveld1993). MEDLL is a parametric mitigation technique that uses the maximum likelihood approach to estimate the amplitude, delays and phases of both LOS and multipath components (Bhuiyan et al., Reference Bhuiyan, Lohan and Renfors2008). In order to determine the shape of the multipath distorted correlation function, the MEDLL uses multiple correlators and a reference correlation function. The MEDLL is an iterative algorithm that subtracts each estimated multipath correlation function from the measured correlation function. After this iterative process is complete, the remainder is the estimate of the LOS correlation function.
There are various implementations of the MEDLL technique (Gao et al., Reference Gao, Liu and Long2013). In this paper, we choose the non-coherent MEDLL for comparison analysis. The non-coherent MEDLL algorithm steps can be summarised as follows (Bhuiyan et al., Reference Bhuiyan, Lohan and Renfors2008):
1. Find the maximum peak of the correlation function from Equation (4) and determine the corresponding delay
${\hat \tau _1}$, amplitude ${\hat \alpha _1}$, and ${\hat \varphi _1}$.2. Subtract a reference correlation function at estimated parameters from the previous step (i.e.
${\hat \alpha _1}\; {{\rm {\cal R}}_{ref}}(t - {\hat \tau _1}){e^{j{{\hat \varphi} _1}}}$ from the measured correlation function ${\hat {\cal R}}$ as follows:
(14)$${{\rm {\cal R}}^1}\left( \tau \right) = \; {\hat {\cal R}} - \; {\hat \alpha _1}\; {{\rm {\cal R}}_{ref}}(t - {\hat \tau _1}){e^{\,j{{\hat \varphi} _1}}}$$3. Determine the new peak of the remaining function
${{\rm {\cal R}}^1}\left( \tau \right)$ at its corresponding delay, amplitude and phase.4. The two previous steps are repeated until the remaining function is below a certain threshold (e.g. set from 0·3 to 0·4 which is used in this paper) or until adding a new delay that does not reduce the mean square error between the measured correlation function and the estimated correlation function.
3. METHODOLOGY
In this paper, the Fast Orthogonal Search (FOS) algorithm (Korenberg, Reference Korenberg1989b; Reference Korenberg1989a; Youssef, Reference Youssef2008) is proposed to model the distorted correlation function to estimate TOA of the LOS and multipath signals. The motivation behind the development of the multipath mitigation technique based on FOS was to find an algorithm that improves the code delay estimation inside the tracking loop of GPS receivers especially for urban canyon navigation environments. The main idea is utilising the high-resolution estimation capabilities of FOS to mitigate the mystical multipath signals. The general architecture of proposed technique inside the delay-tracking loop is shown in Figure 1.
Figure 1. General architecture of FOS-based multipath mitigation technique inside the delay tracking loop.
As shown in Figure 1, after the essential front-end processing and after the carrier has been wiped off, the received signal processing passed into a bank of correlators. The Numerical Control Oscillator (NCO) and PRN generator block produce a group of early and late versions of replica codes. In the case of the classical DLL, the received signal is correlated with each replica in the bank of correlators. Several code tracking techniques (called discriminator in Figure 1) use the output of the correlator bank values as input to generate the estimated LOS delay as output (Bhuiyan and Lohan, Reference Bhuiyan and Lohan2010). A loop filter then smoothes the output of the discriminator.
The proposed algorithm in Figure 1 (highlighted in red) uses the corresponding output of the bank of correlators as input to the FOS algorithm to model the measured correlation function (i.e. distorted correlation function due to multipath effects) and estimate the multipath parameters along with the LOS signal.
The FOS algorithm utilises an arbitrary set of non-orthogonal candidate functions P m(n) and finds a functional expansion of a measured correlation function 
${\hat {\cal R}}(n)$ in order to minimise the Mean Squared Error (MSE) between the measured correlation function and the functional expansion.
The functional expansion of the input 
${\hat {\cal R}}(n)$ in terms of the arbitrary candidate functions P m(n) is given by:
$${\hat {\cal R}}\left( n \right) = \; \mathop \sum \limits_{m = 0}^M {a_m}{P_m}\left( n \right) + e\left( n \right)$$where a m is the set of weights of the functional expansion, P m(n) are the model terms selected from the set of candidate functions and e(n) is the modelling error.
These model terms are a set of GPS L1 ideal correlation functions 
${{\rm {\cal R}}_{ideal}}(n)$ for the middle/prompt correlator of a certain code-delay window range with several multipath delays, phases and amplitudes thereof:
$${P_m}(\alpha, \tau, \; \varphi ) = {\alpha _m}\; {{\rm {\cal R}}_{ideal}}({\tau _m}){e^{\,j{\varphi _m}}}$$where α is the relative amplitude of l-th multipath, τ is the relative delay and φ is the relative phase with respect to LOS signal. These candidate functions are generated offline and saved in a lookup table in memory. The set of multipath parameters can be defined as follows
$$0 \le \alpha \le 1\; $$
$$0 \le \tau \le 2{T_c}$$
$$0 \le \varphi \le 2\pi $$where T c is GPS C/A code chip duration. The code-delay window range is determined based on the estimated correlation peak which, theoretically, can be anywhere within the code delay window range of ± τ w chips (Bhuiyan and Lohan, Reference Bhuiyan and Lohan2010). The code delay window range essentially depends on the number of correlators (i.e. M) and the spacing between the correlators (i.e. Δ) according to
$${\tau _w} = \pm \displaystyle{{(M - 1)} \over 2}\; \Delta $$Choosing non-orthogonal candidate functions does not yield a unique solution for Equation (15). However, FOS may model the measured correlation function with fewer model terms than an orthogonal functional expansion (Korenberg, Reference Korenberg1989a).
FOS begins by creating a functional expansion using orthogonal basis functions such that:
$${\hat {\cal R}}\left( n \right) = \; \mathop \sum \limits_{m = 0}^M {g_m}{w_m}\left( n \right) + e\left( n \right)$$where w m(n) is a set of orthogonal functions derived from the candidate functions P m(n), g m is the weight and e(n) is an error term. The orthogonal functions w m(n) are derived from the candidate functions P m(n) using the Gram–Schmidt (GS) orthogonalisation algorithm. The orthogonal functions w m(n) are implicitly defined by the Gram–Schmidt coefficients α mr and do not need to be computed point-by-point (Korenberg, Reference Korenberg1988; Osman et al., Reference Osman, Nourledin, El-Sheimy, Theriault and Campbell2009).
The Gram–Schmidt coefficients α mr and the orthogonal weights g m can be found recursively using the equations (Korenberg, Reference Korenberg1989a):
$${g_m} = \; \displaystyle{{C(m)} \over {D(m,m)}},\; m = 0,1,\; \ldots, \; M$$where
$$D\left( {0,0} \right) = \; 1$$
$$D\left( {m,0} \right) = \; \overline {{P_m}\left( n \right){P_0}(n)}, \; m = 0,1,\; \ldots, \; M$$
$$D\left( {m,r} \right) = \; \overline {{P_m}(n){P_r}(n)} - \; \mathop \sum \limits_{i = 0}^{r - 1} {\alpha _{ri}}D\left( {m,i} \right),\; m = 0,1, \ldots, M;r = 1,2, \ldots, m$$Further,
$${\alpha _{mr}} = \displaystyle{{\overline {{P_m}(n){w_r}(n)}} \over {\overline {w_r^2} (n)}} = \; \displaystyle{{D(m,r)} \over {D(r,r)}}$$
$$C\left( 0 \right) = \; \overline {{\hat {\cal R}}\left( n \right){P_0}(n)} $$
$$C\left( m \right) = \; \overline {{\hat {\cal R}}\left( n \right){P_m}(n)} - \; \mathop \sum \limits_{r = 0}^{m - 1} {\alpha _{mr}}C(r)$$In its last stage, FOS calculates the weights of the original functional expansion a m (Equation (15)), from the weights of the orthogonal series expansion g m and Gram–Schmidt coefficients α mr (Korenberg, Reference Korenberg1988; Osman et al., Reference Osman, Nourledin, El-Sheimy, Theriault and Campbell2009). The value of a m can be found recursively using
$${a_m} = \; \mathop \sum \limits_{i = m}^M {g_i}{v_i}$$where v m = 1 and
$${v_i} = - \mathop \sum \limits_{r = m}^{i - 1} {\alpha _{ir}}{v_r},\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; i = m + 1,\; m + 2,\; \ldots, \; M$$ FOS requires the calculation of the correlation between the candidate functions and the calculation of the correlation between the input and the candidate functions. The correction between the measured correlation function and the candidate function 
$\overline {{\hat {\cal R}}\left( n \right){p_m}(n)} $ is typically calculated point-by-point once at the start of the algorithm and then stored for later quick retrieval.
The MSE of the orthogonal function expansion has been shown to be
$$\overline {{\varepsilon ^2}\left( n \right)} = \; \overline {{{{\hat {\cal R}}}^2}(n)} - \; \mathop \sum \limits_{m = 0}^M g_m^2 \overline {w_m^2 (n)} $$It then follows that the MSE reduction given by the mth candidate function is given by
$${Q_m} = \; g_m^2 \; \overline {w_m^2 (n)} = \; g_m^2 \; D(m,m)$$The candidate with the greatest value for Q is selected as the model term, but optionally its addition to the model may be subject to its Q value exceeding a threshold level (Korenberg, Reference Korenberg1989a; Osman et al., Reference Osman, Nourledin, El-Sheimy, Theriault and Campbell2009). The residual MSE after the addition of each term can be computed by
$$MS{E_m} = \; MS{E_{m - 1}} - {Q_m}$$The search algorithm may be stopped when an acceptably small residual MSE has been achieved (i.e. a ratio of the MSE over the mean squared value of the measured correlation function (Korenberg and Adeney, Reference Korenberg and Adeney1998) or an acceptably small percentage of the variance of the time-series being modelled). The search may also stop when a certain number of terms have been fitted. In this research, we assumed that the total number of multipath components does not exceed four. Another stopping criterion is when none of the remaining candidates can yield a sufficient MSE reduction value (this criterion would be representative of not having any candidates that would yield a MSE reduction value greater than the addition of a white Gaussian noise series). The order of floating point operations required for FOS is CN + CM 2, where C is the number of candidates searched, n is the number of points in the record and M is the total number of terms constructing the model. The delay of the FOS-based algorithm depends on a number of factors including the data rate, the implementation environment and the hardware used.
In order to estimate the LOS delay, in this research, we assumed that the LOS is present and it has the shortest time of arrival from the satellite to the receiver (shortest transit time). In case the LOS is completely blocked by any obstacle, the proposed algorithm will use the shortest delay estimate. We also assumed that the Doppler shift f D is correctly estimated by the carrier-tracking loop and that all multipath components experience similar Doppler shift.
The performance of the proposed method and the different considered multipath mitigation techniques are compared for two path static channels with path amplitude values 1 and 0·8. The performance comparison of the algorithms mainly focuses on the short-delay multipath scenario because this scenario is generally the most challenging situation. Figure 2 shows the Multipath error envelopes with respect to multipath delay for two static path channels with path amplitude values 1 and 0·8. The upper envelope is obtained for in-phase paths and the lower envelope is obtained for 180° phase shift between the LOS and Non-Line of Sight (NLOS) paths.
Figure 2. Multipath error envelopes with respect to multipath delay for 2 path static channels with path amplitude values 1 and 0·8.
The EML (with 0·2 chip early-late spacing), HRC (with 0·1 chip early-late spacing and 0·2 chip very early-very late chip spacing) and MEDLL are used for the performance comparison. As shown in Figure 2, the proposed algorithm outperforms other methods for the closely spaced multipath scenario. It is also remarkable that the MEDLL algorithm outperforms both EML and HRC multipath mitigation techniques.
4. SIMULATION ENVIRONMENT
In this paper, the SPIRENT GSS8000 simulator, which has the capability to generate realistic simulation multipath environments through an advanced multipath model implemented in the SPIRENT SimGEN™ software, is used to evaluate the performance of the proposed technique under several operation scenarios. Many papers in the literature (Boulton et al., Reference Boulton, Read, Macgougan, Klukas, Cannon and Lachapelle2002) discuss the matching of the simulated multipath environments using the SPIRENT simulator to real multipath environments.
Two different multipath environments are simulated using the SPIRENT GSS8000 simulator using the Land Mobile Multipath (LMM) model (SPIRENT, 2011). The first is a LMM with a rural multipath environment and the second is a LMM with an urban canyon multipath environment. The Rician fading model describes the fading on the LOS signal, however, a Rayleigh model describes the fading on the multipath signals (SPIRENT, 2011). The characteristics of the Rician and Rayleigh models are functions of satellite elevation angle. The details of the Rician and Rayleigh fading models are as follows.
4.1. Rician LOS Fading Model
The Rician model can be defined as (SPIRENT, 2011):
$${\,f_{Rician}} = \left\{ {\matrix{ {2Kv\; exp \lfloor\!\!- K\left( {{v^2} + 1} \right)\rfloor\; {I_0}(2Kv),\; \; x \ge 0} \cr {\hskip10.2pc0,\; \; x \lt 0} \cr}} \right.$$where v is the received signal amplitude relative to the direct path, K is the ratio of direct to multipath power received and I 0 is a 0th order modified Bessel function of the first kind.
4.2. Rayleigh Multipath Model
According to SPIRENT (2011), a modified Rayleigh model describes the fading on multipath channels using three functions. There is a deterministic mean power function, an amplitude noise function following Rayleigh distribution and a delay on multipath channels following exponential distribution. The deterministic mean power reduction in addition to Rayleigh noise can be written as (SPIRENT, 2011):
$${P_h}\left( \tau \right) = \; {P_h}\left( 0 \right) - d\; \times \; \tau \; dB$$where P h(0) is the echo power for zero echo delay (dB), d is the power decay (dB.μs−1) provided by a look-up table and τ is the delay of the echo signal.
The amplitude noise on the multipath channel is randomly computed from a Rayleigh distribution given by:
$${\,f_{Rayleigh}}\left( v \right) = 2\; K\; v\; {e^{ - K{v^2}}}$$where K is the ratio of direct to multipath power received and v is the received voltage relative to the direct path.
The delay on the multipath channel is calculated at random with the exponential distribution given by:
$${\,f_{exp}}\left( \tau \right) = \; \displaystyle{1 \over b}\; {e^{\displaystyle{{ - \tau} \over b}}}$$where b is exponential delay (μs) and it is taken from a look-up table.
The operation of the multipath models is controlled by the contents of a number of Look-up-Tables (LUTs), determined by satellite elevation angle and iteration rate (equivalent to simulated user speed of travel). There is a LUT for every environment type. Tables 1 and 2 show the LUT of rural and urban canyon simulation environments, respectively.
Table 1. Look-up table of rural simulation environment.
Table 2. Look-up table of urban canyon simulation environment.
In order to examine the performance of the proposed FOS-based multipath mitigation technique and compare its performance to the considered multipath mitigation techniques, a test consisting of a set of static and dynamic realistic simulation scenarios was conducted on a SPIRENT GNSS8000 simulator controlled by SimGEN™ software in the TECTERRA Geomatics laboratory at Calgary, Alberta, Canada. Figure 3 shows the experimental hardware setup inside the TECTERRA lab. The reference location was chosen at a start point at Kingston, Ontario, Canada at latitude 44° 13·726′, longitude −76° 2794′ and height 100 m. The experiment duration was 15 minutes.
Figure 3. Experimental setup inside TECTERRA laboratory.
GPS In-phase/Quadrature (I/Q) raw measurements were logged using NovAtel FireHose front-end. The specifications of the digitised signals collected by the NovAtel FireHose front-end are shown in Table 3. The NavINST research group's modified software receiver processed the raw GPS I/Q samples (Borre et al., Reference Borre, Akos, Bertelsen, Rinder and Jensen2007).
Table 3. Settings adopted for data collection.
5. RESULTS AND DISCUSSION
Figure 4 illustrates a frequency domain, time domain and histogram of the raw GPS L1 I/Q samples collected data, respectively. In frequency domain, it clearly shows that the FireHose front-end rotates the received frequencies to the baseband (NovAtel, 2013). In the histogram, it is clear that all four bits of the Analogue to Digital Converter (ADC) are being produced based on the 16 levels present within the histogram.
Figure 4. Frequency domain, time domain and histogram of the raw GPS L1 I/Q samples collected data.
There are nine GPS satellites available above a 5° elevation mask at the initialisation location in the open sky scenario as shown in Figure 5. The open sky scenario is a built-in scenario inside the SPIRENT SimGen™ simulator software. The sky plot shown in Figure 5 illustrates the GPS satellite availability before applying the LMM models.
Figure 5. GPS satellite availability at the initialisation location at open sky scenario.
Two different realistic simulation scenarios including rural and urban canyon were considered. The details about these multipath models are discussed in Section 4.
5.1. Rural Simulation Environment
Figure 6 shows the sky plot of GPS satellites in the realistic rural simulation environment. There are seven GPS satellites available. The low elevation satellites such as PRN 04 and 31 (shown in red in Figure 6) are completely blocked due to the simulated environmental mask. According to the LMM category mask of the rural environment, the signal received from high elevation GPS satellites such as PRNs 05 and 30 are only LOS signals (shown in green in Figure 6). Only the signal received from GPS satellites PRNs 02, 10, 12, 24 and 29 consists of LOS and Multipath components (shown in blue in Figure 6). These signals follow the fading models mentioned in Section 4. The following analysis mainly focuses on the received signals, which is affected by multipath errors. The number of satellites acquired by the GPS software receiver is shown in Figure 7.
Figure 6. Skyplot of GPS satellites at rural simulation environment.
Figure 7. Acquisition results of the GPS software receiver, rural environment.
Figure 8 shows an example of the received correlation function of satellite PRN 10, the estimated correlation function and LOS and NLOS delays of a GPS satellite signal using the proposed method. In order to estimate the delay and phase of the LOS and NLOS signals, Fast Orthogonal Search (FOS) was used to model the measured correlation function by defined pre-designated candidate functions and iteratively adding the term that lowers the MSE of the model by the greatest amount (Korenberg, Reference Korenberg1988). Figure 9 presents the model-fit MSE reduction Qm versus the number of candidate terms of modelling the measured correlation function of GPS PRN 10.
Figure 8. Example of distorted measured correlation function of GPS satellite PRN 10, estimated correlation function and the estimated LOS and NLOS delays using the proposed technique.
Figure 9. Model-fit MSE reduction Qm versus the number of candidate terms of modelling the measured correlation function of GPS satellite PRN 10.
Since the shortest time of arrival from the satellite to the receiver is obtained from the outputs of the proposed method, the pseudorange error of the corresponding PRN is calculated. Figure 10 shows the Root-Mean-Square (RMS) values of the pseudorange estimation errors computed by the proposed method for the satellites PRNs 02, 10, 12, 24 and 29 which consist of LOS and Multipath components (shown in blue in Figure 6). These values are compared to the estimation errors calculated by the considered classical and advanced multipath mitigation techniques discussed in Section 2 such as ELS, HRC and MEDLL.
Figure 10. RMS values of pseudorange estimation errors.
The comparison of the pseudorange estimation errors computed by several techniques is shown in Figure 10. The RMS values of the estimated errors produced by the proposed method are considerably smaller than the ones produced by the other methods for all of the PRNs. This improvement in estimation accuracy is the cost of the additional computation load.
The estimated pseudoranges for all visible satellites were then used in the computation of the positions of the receiver. Figure 11 shows the position errors for the proposed method and other methods. The corresponding RMS values are compared in Figure 12. As was expected from the pseudorange errors, the position solution produced by the proposed method is considerably improved compared to the other considered methods.
Figure 11. East and north position errors, rural environment.
Figure 12. RMS values of position errors, rural environment.
5.2. Urban Canyon Simulation Environment
The proposed method is also examined through another realistic environment simulating urban canyons. Figure 13 shows the sky plot of GPS satellites in the realistic urban canyon simulation environment. There are eight GPS satellites available, the low elevation satellites such as PRN 31 (shown in red in Figure 13) are completely blocked due to the simulated environmental mask. According to the LMM category mask of the urban canyon environment, the signal received from high elevation GPS satellites such as PRNs 05, 12, 29 and 30 are only LOS signal (shown in green in Figure 13). Only the signal received from GPS satellites PRNs 02, 10 and 24 consists of LOS and multipath components (shown in blue in Figure 13) and the signal received from GPS satellite PRN 04 consists only of multipath components (i.e. No LOS component, shown in orange in Figure 13).
Figure 13. Skyplot of GPS satellites in the urban canyon simulation environment.
The comparison of the pseudorange estimation errors computed by several techniques for urban canyon simulation environment is shown in Figure 14. It is clearly shown in Figure 14 that the RMS values of the estimated errors produced by the proposed method are considerably smaller than those produced by the other methods for all of the PRNs. The RMS values of the estimator error produced by the HRC are slightly smaller than the ones produced by ELS for PRNs 02, 10 and 30.
Figure 14. RMS values of pseudorange estimation errors, urban canyon simulation environment.
Figure 15 shows the position errors for the proposed method and the considered methods. The corresponding RMS values are compared in Figure 16. The position solution produced by the proposed method is considerably improved compared to the other methods.
Figure 15. East and north position errors, urban canyon simulation environment.
Figure 16. RMS values of position errors, urban canyon simulation environment.
6. CONCLUSION
In this paper, a new high-resolution code-tracking algorithm has been proposed to mitigate multipath effects for GPS receivers. The new algorithm is based on FOS, which is used to model the distorted multipath correlation function and get better TOA estimation for both LOS and multipath signals. The simulation and data processing results demonstrated the effectiveness of the proposed algorithm in reducing the RMS error values of the estimated pseudorange and positioning solutions compared to other classical and advanced multipath mitigation techniques. It is important to mention that this improvement is obtained at the cost of a large computational load.
ACKNOWLEDGMENTS
The authors would kindly like to acknowledge the technical support of TECTERRA geomatics laboratory (Calgary, AB) in conducting the experimental work using SPIRENT simulator GSS8000. The authors would like to thank NovAtel Inc. for facilitating the use of the front-end Firehose GPS receiver.
