BDS B1I multipath channel statistical model comparison between static and dynamic scenarios in dense urban canyon environment

Global Navigation Satellite System (GNSS) multipath channel models are fundamental and critical for signal simulation and receiver performance evaluation. They also aid the designing of suitable multipath error mitigation algorithms when the properties of multipath channel are available. However, there is insufficient existing research on BeiDou Navigation Satellite System (BDS) signal multipath channel models. In this study, multipath channel statistical models are established on the basis of extensive datasets of the BDS B1I signal. A multipath parameter estimation algorithm is designed to extract information of multipath rays from the intermediate frequency data. The delay, power loss, Doppler fading frequency, and lifetime distribution models for static and dynamic vehicle platforms are established and compared, and the effects of the satellite orbit type and platform speed on the models are analyzed. The results reveal the detailed distribution and variation characteristics of the multipath parameters and are valuable for the development of accurate urban navigation systems.


Introduction
Multipath is a major factor that significantly impedes the performance of GNSS receiver. Its reception causes a synchronization bias with respect to the direct Lineof-Sight (LOS) signal for receiver tracking loops, thus leading to errors in code and carrier phase pseudorange measurements. The multipath signal is usually modeled as a delayed and attenuated replica of the LOS signal. Proper understanding of its characteristics in different environments becomes important for the GNSS receiver design and channel simulation.
Conventionally, Code-Minus-Carrier (CMC) observables (Misra and Enge 2006) have been used to isolate multipath code errors and study their characteristics. For a receiver on a static platform, multipath errors exhibit site-specific and periodically repeatable patterns (Axelrad et al. 2005;Wang et al. 2015). However, the CMC method cannot indicate the number of multipath rays and their parameters, such as the delay, attenuation, Doppler frequency, and lifetime. The ray-tracing technique is employed to calculate all possible multipath signals precisely for a specific scenario (Weiss et al. 2007;Gowdayyanadoddi et al. 2015;Lau and Cross 2017;Panicciari et al. 2018). This technique is more applicable to multipath analysis for small-scale scenarios. It encounters difficulties when applied to large-scale environments, such as urban environments, where the precise three-dimensional maps as well as the surface material information are not available or are too complex to assemble.
For large-scale environments, multipath statistical channel models must be utilized. The German Aerospace Center (Deutsches Zentrum für Luft-und Raumfahrt;DLR) conducted a series of research campaigns on the Global Positioning System (GPS) and Galileo navigation satellite system (Galileo) L1

Open Access
Satellite Navigation https://satellite-navigation.springeropen.com/ *Correspondence: xin.chen@sjtu.edu.cn Shanghai Key Laboratory of Navigation and Location-based Service, School of Electronic Information and Electrical Engineering, Shanghai Jiao Tong University, Dongchuan Road 800, Shanghai 200240, China multipath channel, for which the researchers used a Zeppelin to transmit a 100-MHz broadband signal in the L band. The multipath characteristics, such as the delay spread, power attenuation, lifespan, and Doppler frequency, were reported and compared for urban and suburban environments and for land mobile and pedestrian platforms Lehner 2004, 2008;Lehner and Steingass, 2005;Steingass et al. 2017). Falco et al. (2018) implemented this model by a hardwarein-the-loop approach to evaluate the multipath effects in railway environments. However, the consistency of these characteristics with the real GNSS multipath channel still needs to be validated. Moreover, it is not clear whether these models are applicable for eastern Asian metropolitan cities, such as Shanghai, Hong Kong, or Tokyo, because the urban building structures of these cities are quite different from those of Munich.
As an initial effort to study the multipath channel, particularly for BDS B1I signal in eastern Asian cities, a campaign was conducted to establish the channel models using real GNSS signal measurements. A broadband Intermediate Frequency (IF) data recording device was used to collect BDS B1I signals. A Real-Time Kinematic (RTK)-GNSS/Inertial Measurement Unit (IMU) integrated device was employed to obtain the ground truth of the platform positions during the IF data collection. In the first phase of the research campaign, IF data collections of the B1I signal for static platforms in the Shanghai-Lujiazui dense urban canyon area were performed, and the corresponding models for the delay, average power loss, and Doppler fading frequency were preliminarily established (Chen 2018;Wang et al. 2018). In the second phase, the IF data for dynamic platforms in the same area were extensively collected, and the B1I signal multipath models were developed, as discussed in this report.
A major aspect posing a challenge to research on the dynamic platform multipath channel is the highly dynamic channel conditions, which make it difficult to estimate the multipath signal parameters correctly if the conventional multipath estimation methods are adopted, such as the Multipath Estimating Delay Lock Loop (MEDLL) method (Van Nee 1992). Therefore, a multi-parameter estimation method for dynamic channels was carefully designed, including the number of multipath signals and their multipath delays, signal amplitudes, Doppler fading frequencies, multipath lifetimes, and Non-Line-of-Sight (NLOS) signal detection. Finally, statistical channel models of the B1I signal for both static and dynamic platforms were fully established, and the model differences for different BDS constellations and different motion speeds were compared.
These results are useful for simulating high-fidelity urban channels for advanced navigation algorithm developments.
The remainder of this paper is organized as follows. In Sect. 2, the multi-parameter estimation method used for multipath data extraction from the real B1I signal is described. In Sect. 3, the IF data collection system is introduced, and two data collection campaigns performed in the Shanghai-Lujiazui area are described. Then, the BDS B1I signal multipath channel models are established and analyzed in Sect. 4. Finally, it is concluded with a brief summary and a description of planned future work.

Multipath multi-parameter estimation method for dynamic channel
For GNSS radio channel estimation, the multipath signal with a delay within one and a half code chips delay relative to the LOS is more interesting because the multipath signal falling in this delay affects the pseudorange measurement accuracy most significantly. The MEDLL method (Van Nee 1992) with the maximum-likelihood principle was proposed to estimate the delay, amplitude, and phase bias. However, the MEDLL method assumes a stable channel condition and needs three-dimensional exhaustive search, which prohibits its extensive application. The least-squares-based iterative multipath superresolution algorithm (Nam and Kong 2013) employs an iterative update formula to estimate the delay and complex amplitude simultaneously to reduce the computational complexity. However, the update step size must be selected carefully for convergence. The fast-iterative maximum-likelihood algorithm (Sahmoudi and Amin 2008) employs another iterative closed-form solution to estimate multipath delays. The Coupled Amplitude Delay Lock Loop (CADLL) algorithm (Chen et al. 2011(Chen et al. , 2013 utilizes several units of the Delay Lock Loop (DLL) and Amplitude Lock Loop (ALL) to estimate and lock to the multipath delay and complex amplitudes. The coupled DLL and ALL method improves the estimation resolution and tracks multipath parameters sequentially.
However, it is found that there has not been a readily suitable method for GNSS multipath parameter estimation in urban dynamic environments. Thus, a systematic method for estimating multipath delay, amplitude, Doppler fading frequency, lifetime, and the NLOS signal is introduced in this report.
The multipath channel for the GNSS land vehicle application in an urban environment is time variant, and it can be modeled as a superposition of discrete paths, given by (Rappaport 2009). where α l (t) and τ l (t) denote the complex weight and the propagation relative delay of the l-th path, respectively, l = 0 denotes the LOS signal, c(t) is the spreading code of the GNSS signal, f is the carrier frequency of the LOS signal, the Doppler fading frequency v l (t) is the carrier frequency difference of the l-th multipath relative to LOS in which it is caused by the relative movement between the satellite, reflector, and receiver, T δ is the time interval of interest, L is the number of paths, and n(t) is the white Gaussian noise with power spectra equal to N 0 /2 . Without loss of generality, the delays { τ 0 , . . . , τ L−1 } are assumed to be distinct, where τ 0 < τ 1 < · · · < τ L−1 .
The multipath lifetime ε(t) is used to describe the existing time of a multipath ray in a real environment. Therefore, the vector θ l [τ l , α l , v l , ε l ] is defined to characterize a multipath signal uniquely. Then, θ L−1 l=0 θ l denotes the entire multipath parameter space, which is of interest for estimation.
In a GNSS receiver, x(t) in (1) is first down-converted to the IF band, following which it is sampled and correlated with the local replica signal. Let the correlation interval be defined by T = NT s , where T s is the sampling period and N is the total sample number in the correlation interval. The resultant correlation value has the form where R δτ k − τ l,k is the normalized autocorrelation function of the spreading code sequence, n k is the complex white Gaussian noise with zero mean and variance of N 0 /(2T ) , δτ k and δθ k are, respectively, the time and phase synchronization error between the LOS signal and the local replica signal, and τ l,k and ϕ l,k are, respectively, the average relative delay and carrier phase bias of the l-th multipath path with respect to the LOS in the k-th interval for which τ l,0 = 0 and ϕ l,0 = 0.
When T is small, it is reasonable to assume that the channel is invariant in the small interval. Thus, one can have τ l,k ≈ τ l,k , v l,k ≈ v l,k , and α l,k ≈ α l,k . To simplify the notation, the subscript k is omitted in the following expressions; however, it should be noted that the estimated parameters from the correlation value z are considered as in the current interval. To increase the estimation observability, a bank of correlation values with equal spacing intervals in the delay dimension are computed in the receiver, given by (1) where z(k) is the correlation bank vector and In (4), ∆ is the spacing interval between two adjacent correlators, and M is the total number of correlators in the bank.
In this study, on the basis of the CADLL architecture (Chen et al. 2011(Chen et al. , 2013, an algorithm for functioning in a dynamic channel condition was developed to estimate the multipath delay, complex amplitude, and Doppler fading frequency from the raw IF GNSS signal. First, the correlation value vector z(k) in (3) is decomposed into L separate single path correlation value vectors, i.e., where θ p (k − 1) is the estimated parameter vector for the p-th multipath from the previous k − 1 interval and ẑ p k|θ p (k − 1) is the predicted correlation bank vector for the p-th multipath according to the signal model in (2) with the aid of θ p (k − 1): In (6), α p,k−1 , ϕ p,k−1 , and τ p,k−1 are the estimated multipath amplitude, carrier phase bias, and delay of the p-th multipath from the k − 1 interval. The initial multipath parameters can be obtained via a coarse grid search. For each decomposed single path, a narrowspacing DLL together with two ALLs, called a "unit, " are used to estimate and track the delay, the in-phase amplitude, and the quadrature-phase amplitude. The spacing of the DLL is set as 0.1 chip to guarantee a good tracking accuracy for the delay.
Because the delay and amplitude of a multipath signal usually have fast fading variations in a dynamic channel, the original CADLL architecture was modified by incorporating a Kalman filter to track the variations.
The state vector of the Kalman filter is defined as s = τ ,τ , α I ,α I , α Q ,α Q , ϕ, v T , where τ is the derivative of the delay τ , α I and α Q denote the in-phase and quadrature phase amplitudes, respectively, and α I and α Q denote the derivatives of α I and α Q , respectively. The state transition can be derived as where Φ is the transition matrix, A 2x2 = 1 T k 0 1 , and w 8x1 is the state transition noise.

The measurement equation of the Kalman filter is
where y k is the measurement vector, in which τ k is the multipath delay being tracked by the unit's DLL for the lth multipath ray, α I|P k and α Q|P k are the in-phase and quadrature phase prompt correlator values of the unit for the l-th multipath, and ϕ k is the multipath phase bias measurement that is given as arg α Here, H is the observation matrix in which B 1x2 = 1 0 , and r 4x1 is the measurement noise. It is necessary to note that w 8x1 and r 4x1 are affected by the vehicle motion, signal carrierto-noise ratio, and environment type. They can be carefully tuned and stored in a lookup table for future use. As soon as a multipath signal is detected, a unit containing the Kalman filter shown in (7) and (8) is allocated for tracking it. The C/N 0 of this multipath signal, estimated by the narrow wide band correlation method (Parkinson and Spilker 1996), is used as the metric for detecting a loss of lock. A loss of lock occurs for this multipath when its estimated C/N 0 is lower than a threshold. In this case, this multipath is regarded as being vanished. Accordingly, the emerging and vanishing time of the multipath is recorded by the unit, and its lifetime ε can be readily obtained.
Next, the model order determination can be applied to determine whether the input signal is over-fitted, i.e., whether the supposed number of multipath signals is greater than the actual number. Classical information theory methods for model selection can be used to estimate L (Wax and Kailath 1985;Was and Ziskind 1989). According to the method of Wax and Kailath, the optimal model order L is the one that achieves the Minimum Description Length (MDL) given the observed data z(k): where In (10), P Z|θ (L) is the probabilistic model of the observed dataset Z = {z(1), · · · , z(N )} , where N is the number of epochs that are lumped to make the set, θ (L) is the optimal parameter vector estimate under the assumption of L, and χ is the number of free parameters in the vector θ (L) . To compute the MDL value, sample-covariance estimation and eigenvector decomposition are applied. More-detailed information about the MDL method can be found in previous research (Wax and Kailath 1985;Wax and Ziskind 1989). Then, all the valid parameters for each path along with the optimal model order L are grouped together as a new set of multipath estimates and fed to the next epoch.
In this study, a correlator bank with M = 127 is used for each satellite channel. The correlation interval T was set to 10 ms, and the minimal detected multipath power is −150 dBm, which corresponds to approximately maximal 20 dB power attenuation relative to the LOS signal. This minimal multipath power detection level is carefully designed due to the cross-correlation mitigation level between different B1I PRN codes. Multipath signals lower than this power level cannot be guaranteed to enable reliable detection and are therefore treated as noise. The maximally detected multipath number for each satellite is 4, and the maximal delay estimation window is 600 m. The estimation resolution for two closely spaced multipath signals in delay is 0.1 chip, which corresponds to approximately 15 m for the B1I signal.
In the algorithm mentioned above, it is supposed by default that the first arrival path is the LOS signal and all the delays and carrier phase bias are estimated relative to the first arrival path. In the NLOS case, it encounters a strong multipath while the LOS signal is completely blocked by some obstacle. Thus, all the estimated parameters are erroneous because their delays and carrier phases are relative to the multipath signal, which is being tracked by the receiver. Therefore, finding a way to detect and remove the NLOS situation is needed. During the data collection, a precise RTK-GNSS/IMU receiver, e.g., the NovAtel FANS system, was used to work simultaneously and share the same antenna with the signal recording system. Precise trajectory position, vehicle velocity, and satellite ephemeris can be logged, so that a precise predicted pseudorange for each satellite at each epoch is available. Similarly, the pseudorange for the first arrival path in (2) can be measured as well. Regarding these two pseudoranges, if their difference is larger than a threshold, a NLOS case can arise, that is where ρ meas (k) denotes the pseudorange measurement corrected by the ionosphere delay, troposphere delay, satellite clock error, and local clock error for the first arrival path in the current epoch, ρ pred (k) stands for the predicted pseudorange computed with precise receiver position and satellite ephemeris, and the bound ǫ ρ is empirically set as 15 m according to the extensive recorded data. According to the NLOS detection method given in (11), the local clock error must be carefully estimated because it directly influences the accuracy of the measured pseudorange ρ meas (k) . For this, an open-sky environment is used at the beginning of the data collection and continued for a while to let the receiver correctly compute its initial position and the local clock error. Following this, the van is driven to a more densely urbanized environment to collect the multipath data. In the dense urban environment, the uncontaminated LOS signals are selected to maintain local clock error estimation. When no clean LOS signals are available in some particularly poor environments, the local clock error could be treated as constant until the van is driven out of the area. When new, clean LOS signals are observed, the receiver local clock error estimation is continued. It has been proved by experiments that this method can work well in the data collection campaign.
When a NLOS case is detected, the signal tracked by the receiver is just the NLOS signal. Thus, the pseudorange error computed in (11), the amplitudes and Doppler frequency that are estimated in the tracking loop, and the time period during which the NLOS case exists constitute the parameter vector θ l of the NLOS signal, which are the delay, Doppler fading frequency, amplitude, and lifetime. Figure 1 shows the diagram of the multipath parameter estimation method.
It needs to be pointed out that the estimation method proposed is a universal multipath parameter estimation algorithm for dynamic channel conditions. It can be readily extended to other BDS satellite signals, such as B1C, B2a, and B3I, and other GNSS satellite signals.

BDS B1I raw signal collection campaign
Urban canyon area is the main research environment in this study because it is one of the most important scenarios for navigation applications but suffers badly from the multipath phenomenon. The Shanghai-Lujiazui area is chosen as a representative dense urban canyon environment. This district covers approximately 1.7 km 2 , and most of the buildings are from 80 to 200 m high. The main street width in this area is approximately 60 m, and the side street width is 20-30 m. Two data collection campaigns were performed from 2015 to 2018, one of which was for a static platform scenario and the other for a dynamic platform scenario. A data collection system is developed for collecting the real BDS B1I multipath signal in the urban canyon environment. A broadband IF data record device is designed in which the bandwidth is more than 40 MHz and the complex sampling frequency is 62 MHz, covering GPS L1CA and BDS B1I signals. A one-stage direct down-convert in-phase/quadrature orthogonal sampling technique is used in the device design with a local oscillator frequency of 1568 MHz. Thus, the resultant IFs of GPS L1CA and BDS B1I are 7.42 and − 6.902 MHz, respectively. The sample quantization is Signal decomposition U(L-1) DLL DLL ALLi Kalman filter: Kalman filter: ALLi ALLq ALLq 8 bit. This signal recording device is mounted on a van with the antenna fixed on the vehicle roof. The NovAtel FANS system is mounted on the vehicle as well and shares the same antenna with the recording device. This system provides a precise trajectory position log, which is used in the NLOS detection process. The data collection system is shown in Fig. 2.

Campaign I: Data collection on static platform
From 2015 to 2017, a series of real BDS B1I IF raw signal collections were performed for the static platform scenario. Sixteen equally spaced grids were designed in the Shanghai-Lujiazui area, each of which covered approximately 0.1 km 2 . Twenty-six collection spots were employed during this campaign. The distribution of these spots around this area and some land view examples are shown in Fig. 3. In each collection spot, no less than 30 min of raw IF BDS B1I data were recorded.

Campaign II: Data collection on dynamic platform
In 2018, the research group performed extensive onfield data collection in the Shanghai-Lujiazui area for the dynamic platform scenario, as shown in Fig. 4. Approximately 37.7 km of travel distance was covered in this campaign. The vehicle speed varied between 0 and 50 km/h according to the traffic circumstances. The zero speed and the 50 km/h speed correspond to the stop in front of a red traffic light and the highest allowed speed in this area, respectively.
All the raw IF data were processed by the proposed multipath parameter estimation method, and the parameters of these detected multipath signals were recorded and used to establish the model.  multipath delays and corresponding C/N 0 over time. The x-axis is time (s), the y-axis is multipath delay (m), and the C/N 0 is shown with colors. The redder the color, the higher the C/N 0 . The strip at delay 0 in this chart is the LOS signal, whose C/N 0 stays between 42 and 45 dB·Hz. The C/N 0 of most multipath signals varies between 30 and 35 dB·Hz. It is noticed that four multipath segments can be observed in the chart. Their average delays are approximately 100, 110, 135, and 15 m, respectively, for multipath segments 1 through 4. The fluctuations for both the C/N 0 and delay around their average are caused by noise. For an ideal static multipath reflection environment with a GEO satellite signal, it is commonly believed that the multipath signal basically remains stable. However, the real data show that the multipath signal may have noticeable changes, even within a short time (less than 5 min in this case). This probably is caused by the nearby moving reflectors, such as cars and pedestrians, and the nonhomogeneous materials on the reflection plane. Therefore, one can immediately judge that the multipath delay and C/N 0 evolve with time, but their evolution behavior depends on the specific environment and the satellite constellation. Limited multipath lifetime does exist, even for a GEO signal on a static platform. Figure 6 shows the detected multipath delay and C/N 0 evolutions in a short trajectory data collected in Campaign II for the dynamic scenario. It can be observed that many delays are either increasing or decreasing since the vehicle is moving. These multipath have varying lifetimes, depending on the geometric relationships between satellite, receiver, and reflectors.

BDS B1I signal multipath channel statistical model establishment and analysis
In this section, the multipath statistical models are derived and analysed based on the real multipath data obtained from the two collection campaigns. Since the vector θ [τ , α, v, ε] uniquely characterizes the behavior of a multipath signal, the focus is on the modeling of the four core multipath parameters: the delay distribution for τ , the power loss delay profile for α , the Doppler fading frequency distribution for v , and the lifetime distribution for ε . The multipath data from two collection campaigns, the static platform scenario and the land vehicle dynamic platform scenario, are used for establishing the corresponding models. The effects of the satellite constellations and the vehicle speeds on these models are analyzed.

Multipath delay distribution model
The multipath delay distribution shows the probability that a multipath appears with a certain delay with respect to the LOS in a specific environment (for example, an urban canyon). This model provides insight into how multipath signals spread along the delay dimension if simulating a multipath scenario.  7 and 8 illustrate the multipath delay histograms of the B1I signal from three different satellite orbit types-Medium Earth Orbit (MEO), Inclined Geostationary Satellite Orbit (IGSO), and GEO-that are extracted from the static platform Campaign I data and the dynamic platform Campaign II data, respectively. To build these histograms, all possible multipath delays from 10 m (the minimum) to 600 m (the maximum) are partitioned into evenly spaced small intervals, each of which is 5 m long. Then, the number of multipath whose delays fall into each interval is counted. All the histograms except for the GEO multipath in the static scenario in Fig. 7 have a similar distribution trend that can be fitted by utilizing a gamma function. The distribution shapes do not differ much between different orbit types or movement status. The dissimilar histogram for the GEO multipath in the static scenario is caused by the small number of multipath extracted from the campaign data. If the average delay, 141 m, is compared with the average delay, 155 m, of the GEO multipath collected from the dynamic scenario campaign data, as shown in Fig. 8, one finds that they are close to each other. Therefore, this result indicates that the delay distribution depends mainly on the surrounding environment of the receiver.
Furthermore, the effects of the vehicle speed on the delay distribution are analyzed. Six groups of multipath delays of the dynamic collection data are extracted according to the vehicle speed: 1-5, 5-10, 10-15, 15-20, 20-25, and 25-30 km/h. The delay distribution histograms are displayed in Fig. 9. A quick sanity check shows that all these six histograms have a similar distribution, and their average delays τ s are close to each other. This result verifies again that the delay distribution is determined by mainly the environment other than the receiver moving speed.
Thus, the multipath delay distribution, as an examination of Figs. 7, 8 and 9, is given by where a τ and b τ are the model parameters, and Γ (a τ ) is the gamma function.

Multipath power loss model
The multipath power loss model shows how the mean multipath power varies with respect to multipath delay. The bistatic radar model can be used to predict the received multipath power because of the scattering in the far field (Seidel et al. 1991;Van Rees 1987). The bistatic radar equation, which describes the radio wave traveling in free space, impinging on a distant scatter object and then reradiating in the direction of the receiver, is given by where P i is the power of the incident wave at the reflection/scattering object, P mp is the multipath power received by the receiver, and the Radar Cross Section (RCS) is defined as the ratio of the power density of the signal reflected/scattered in the direction of the receiver to the power density of the radio wave incident upon the (13) P mp = P i + RCS − 20nlog(τ ), object. The RCS value accounts for the reflection/scattering loss at the object surface of two different media. The term 20nlog(τ ) indicates the propagation loss from the reflection object to the receiver, where n is the decay coefficient. It can be predicted from (13) that near large objects usually produce a stronger multipath than small far objects due to the different RCS values and the excessive traveling path τ.
Hereinafter, the focus is on the multipath attenuation, 20log α (dB), with respect to the LOS instead of the multipath power to isolate the effects of the LOS power. Based on the statistical delay analysis of Figs. 7, 8 and 9, the multipath power attenuations in each bin are averaged to obtain the power loss model plot. Figure 10 shows the multipath average power attenuation versus the delay for static platform Campaign I data (left) and dynamic platform Campaign II data (right). It is observed that the multipath average power attenuation versus the delay has an approximately linear decreasing trend whether the receiver is static or dynamic. The average power attenuation PL(τ ) can be well approximated by the form where PL 0|dB = − 12.35 dB and d dB = − 0.0016 dB/m are the model parameters obtained from the fitted curve. Comparing (13) and (14), one finds that they are basically consistent. Both models have a constant term, RCS in (13) and PL 0|dB in (14), and they are irrelevant with the delay τ . In addition, the power loss for both models has a decreasing trend as the delay τ increases. However, the power loss PL(τ ) in (14) follows a linear decline with respect to τ , while the one in (13) declines linearly with respect to log(τ ). This difference is because the model in (13) is usually used in long-distance scattering power loss prediction, e.g., a few kilometers or tens of kilometers. The model in (14) is applied to a short distance, e.g., 0-600 m. In this small distance, the average multipath power loss is almost the same but with a slightly declining trend.

Multipath Doppler fading frequency model
The multipath Doppler fading frequency is generated by the relative geometry change between the satellite, reflector, and receiver. It affects the periodic behavior of observation errors caused by multipath. For convenience, the term "fading frequency" or "fading" is used to denote the multipath Doppler fading frequency in the rest of this paper. Figure 11 shows the multipath fading frequency histograms for B1I MEO, IGSO, and GEO signals obtained from the static platform Campaign I data. It is observed that all the fading frequency histograms are approximately symmetric around zero, and their absolute values are mainly less than 1 Hz. Their distribution can well be fitted by an absolute exponential function, which is given by where P s fading (v) stands for the distribution probability of the fading frequency for the static platform scenario, and µ is the model parameter. It is found that the fading frequency of the GEO multipath is very close to zero because its satellite movement is almost stable. Figure 12 shows the multipath fading frequency histograms for B1I MEO, IGSO, and GEO signals obtained from the dynamic platform Campaign II data. The Gaussian function fitting curve and the absolute exponential function fitting curve are also illustrated for each plot. The multipath fading frequencies in the dynamic scenario can be up to tens of hertz, and they are better fitted by a zero mean Gaussian function where P d fading (v) stands for the fading frequency distribution probability for the dynamic scenario, and σ v is the standard deviation of the Gaussian function. The distribution difference of the multipath fading frequency for the two scenarios occurs because, in a static or quasistatic scenario, a small fading frequency is generated by the slow relative geometry change caused by the satellite. Therefore, most frequencies stay around 0 Hz, thus leading to an exponential-shape-like distribution. When the receiver has substantial movement, the multipath fading frequency starts spreading toward two sides, and a Gaussian-like distribution gradually becomes a superior model.

Multipath lifetime model
The multipath lifetime parameter, denoted by lft or ε , reflects the time duration for a multipath from emergence to extinction. It is useful for simulating a varying multipath scenario. The factors that affect the lifetime of a multipath ray are threefold: the limited size of the reflecting object, the relative movement between the satellite, reflection source, and antenna, and the heterogeneous materials of the same reflection plane.
Because the multipath parameter estimation method can detect and track a multipath ray, the time duration of a multipath is able to be measured. The lifetime measurements of all multipath are sorted into three groups according their orbit type: the MEO, IGSO, and GEO groups. For each group, the distribution is obtained from the histogram with a bin size of 1 s. Figure 13 shows the probability histograms and the Cumulative Density Function (CDF) curves of the three groups for the static scenario Campaign I data. It is found that the multipath   (Chen 2018) signals of the three orbit types have a large portion of lifetime less than 10 s, although the experiment is on a static platform. The large number of short lifetime multipath is probably caused by the reflection generated from moving objects in the busy traffic in this scenario. Furthermore, the Probability Density Function (PDF) histogram of the lifetime is more spread toward a long tail because the variation of the satellite movement is smaller. Finally, the multipath lifetime of the GEO signal is significantly longer than those of the MEO and IGSO signals. Figure 14 shows the multipath lifetime histograms and CDF functions for B1I MEO, IGSO, and GEO signals from the dynamic Campaign II data. Three values are used to reflect the distribution characteristics: α(95%) is the value for which cdf(α) > 95%, α(50%) is the value for which cdf(α ) > 50%, and the probabilistic average value of lifetime lft is defined as lft = N n=1 lft n p n , where lft n is the lifetime value at the n-th bin, and p n is the possibility of the corresponding bin. Two candidate probability functions are used to model the lifetime distribution: the exponential function (green line) and the Generalized Pareto Distribution (GPD) function (pink line). Both functions have been used to model lifetime distributions in statistical analysis. The exponential function is a special case of the GPD function. Figure 14 shows that the lifetime distributions of all three orbit satellite signals on a dynamic platform have very similar shapes and magnitude levels. This results obviously differ from those of the static platform data shown in Fig. 13, which means that platform movement is the major factor that determines the lifetime distribution instead of satellite movement. Moreover, the GPD function fits the histogram better than the exponential function.
Because the orbit type does not affect the lifetime distribution for the dynamic scenario, all multipath data of the dynamic Campaign II data are separated into six groups (0-3, 3-10, 10-15, 15-20, 20-25, and 25-30 km/h) according to their vehicle speed. Next, the statistical values α(95%) , α(50%) , and lft are computed for each group. The effects of different speeds on the lifetime distribution can thus be analyzed. The results are shown in Fig. 15. One can observe that multipath lifetime distribution remains approximately stable when the velocity exceeds 10 km/h. Finally, the analytic multipath lifetime distribution is derived by using the GPD function according to the results shown in Figs. 13, 14 and 15.
(17) where ξ ε is the shape parameter, σ ε is the scale parameter, and η ε is the location parameter. All the multipath parameter models, along with their corresponding model parameters at different scenarios, orbit types, or vehicle speeds, are summarized in Table 1 of the Appendix.

Conclusion
In this study, the real BDS B1I data collection campaigns and the corresponding multipath channel modeling for urban environment are presented. Based on these results, the statistical models of the multipath delay, power-delay loss, Doppler fading frequency, and lifetime are analyzed and compared. A summary of the major findings in this study is listed below.
(1) The gamma function fits the distribution of the multipath delays. The satellite elevation other than the satellite orbit type or the vehicle speed affects the average delay and the delay model parameters.
(2) The multipath signal power has a linear average power-delay model that is not affected by the satellite orbit type, vehicle speed, or satellite elevation. (3) An absolute exponential distribution fits the fading frequency distribution of the multipath signal for a static scenario. In this case, the satellite orbit type affects the model parameters significantly. The GEO multipath has an order of 0.001 Hz fading frequency, while the MEO multipath has an order of 0.1 Hz. The IGSO multipath fading frequency lies between them. For a dynamic scenario, the Gaussian model fits the Doppler fading frequency distribution better, and no obvious difference can be seen for the satellite orbit types. (4) The GPD function fits the distribution of the multipath lifetime. The α(95%) multipath lifetime can be up to 720, 87, and 26 s for GEO, IGSO, and MEO, respectively, in a static scenario. However, the α(95%) multipath lifetimes of the three orbit types do not have a significant difference in a dynamic scenario. The vehicle speed affects the lifetime distribution parameters.
These findings provide insights into the multipath channel characteristics. In the future, it is planned to expand the research to other BDS satellite signals, other constellation systems, and other typical environments, such as tree-lined streets or the lanes down elevated roads.  Table 1.