Abstract
A novel MIMO (multiple input multiple output) satellite channel model that allows the generation of associated channel impulse response (CIR) time series depending on the movement profile of a land mobile terminal is presented in this paper. Based on high precise wideband measurements in L-band the model reproduces the correlated shadowing and multipath conditions in rich detail. The model includes time and space variant echo signals appearing and disappearing in dependence on the receive antenna position and movement, and the actual azimuths and elevations to the various signal sources. Attenuation and path delays relative to the hypothetical line of sight (LOS) ensure usability for ranging purposes. Parameters for car and pedestrian applications in urban and suburban environments are provided. The channel characteristics are determined independently of the transmitted signal. Therefore the usability, for example, for GPS and GALILEO, as well as wideband communication services from hovering platforms, is given.
1. Introduction
Augmentation systems are widely used to improve the localization accuracy of global navigation satellite system (GNSS) receivers by providing estimates of most error sources. But one source of error remains challenging in land mobile applications: shadowing and multipath which are caused by the specific surrounding of a receiver. Depending on delay, phase, and power of each echo signal, the autocorrelation function, which is evaluated for the purpose of ranging, is distorted. The resulting ranging error can largely vary in such an environment [1], even when the receiver moves by just a few meters.
Due to the low receive power of satellite signals it is particularly difficult to determine the precise channel state. New promising receiver architectures are based on estimating the AOA (angle of arrival) of echo signals, for example, with antenna arrays [2], or estimating and tracking the phases and delays of echoes with, for example, particle filters [3]. For the later one a movement model and the geometric representation of reflectors are required. The estimation performance can be further improved by taking correlations among the different satellite links into account. In order to develop and validate such new receivers, wideband channel models are required that implement the urban electromagnetic wave propagation behavior in adequate detail.
Such models are as well required in space/air-to-ground communications, where the most limiting factor is the presence of obstacles blocking the LOS. One of the few alternatives to overcome this deleterious effect is transmitter diversity, which improves link availability. Its performance is directly related to the degree of spatial correlation between the alternative signal and echo paths. Moreover, there is growing demand in adding more precise ranging capability to communication systems, for example, to support urban personal navigation. One example is the use of synchronization and pilot symbols in LTE [4], which is operated in L-band in several countries. In scenarios with nearby base stations or links from/to hovering platforms propagation conditions are within the focus of this work.
Realistic regeneration of space-time variance of echo signals and execution time are key parameters for simulations with wideband models. Traditional approaches are either purely statistical or deterministic. Statistical models recently presented in [5–7], which extend the original concept of [8], are very fast with low complexity; however the delay and phase evolution over time and space of single echoes and the AOA information cannot be deduced. On the other hand ray tracing or physical optics methods like in [9] are very accurate and provide all the necessary details for the ranging error problematic in GNSS. But the generation of complex environments with house facades including, for example, window sills, frames or rain pipes, streets with vegetation, traffic signs, poles, wires, and potentially moving objects like cars and bikes, are exhaustive to create and the processing time grows enormously for such a level of environmental detail.
Hybrid models are fast and accurate: they combine a virtual city with simplified electromagnetic interactions and statistical distributions. Such approaches do exist for ground-to-ground communications, for example, the WINNER II [10] channel model, which is a stochastic model that even takes clusters of reflections into account. One drawback of this model is that the reflector distribution is not validated by measurements. Moreover it is limited to “terrestrial” transmit antenna heights. In case of air/space-to-earth propagation a narrow band measurement-based Markov model was presented in [11], but this model does not provide delay or amplitude information of single echo signals. Some temporal and spatial correlation properties from S-band measurements are discussed in [12]; however these results were obtained at a specific elevation of about 20° and do not allow modeling the strong elevation dependency reported in [13]. Also important to note here is the fact that none of these models considers the excess delay of the first propagation path in non-LOS condition, which is one of the most essential parameters to be accounted for in ranging based positioning. The physical-statistical MIMO channel for LMS communications presented in [14] is based on a cluster concept and uses the same methodology as in [15, 16]. The major drawbacks of these models are vague and very simplified scatter distributions and backscatter properties and the fact that variations in the environment, for example, movements, are not accounted for.
The goal of the wideband MIMO satellite channel model presented in this paper is to take high detailed environmental and geometry related deterministic effects, as well as correlations among the different satellite signal paths into account, in order to develop more effective multipath mitigation techniques and to increase service continuity and quality for mobile receivers in urban and suburban environments in the future. In absence of high delay resolution satellite channel measurements in the past, we performed a unique channel sounding campaign where we used a Zeppelin to simulate a satellite, outlined in Section 2. In Sections 3 and 4 a channel parameter estimation is applied and the relevant effects and statistics are analyzed. Finally the resulting channel model that combines statistically relevant data with observed deterministic effects is described and the generated output is exemplified in Sections 5 and 6, respectively.
2. Land Mobile Satellite Channel Measurement Campaign
The main measurement requirements were realistic far field propagation conditions, a delay resolution in the order of a few ns, and a dynamic range of at least dB to detect all the relevant echoes that can affect high sensitive receivers today. To fulfill these requirements a wideband channel sounder transmitter TX on a flying platform was selected for the setup. Positioned in a distance of about 5 kilometers from the channel sounder receiver RX, which was located in the multipath environment, the received power was large enough to detect the relevant echoes and the wavefront to the zone of reflecting objects around the receiver was almost parallel, as the radius of this zone is typically about m. To precisely measure the propagation effects, TX and RX were synchronized by use of customized rubidium atomic clocks that provided an Allan standard deviation below s even when exposed to vibrations and fast temperature changes on board of the Zeppelin’s cabin.
2.1. Measurement Setup
The sounding signal was transmitted from the hovering Zeppelin at MHz with 10 W equivalent isotropic radiated power and a bandwidth of MHz (see Figure 1). The RX data was collected on the measurement van, which also carried two video cameras, one with front view on the street and the second one with a fish eye lens on the roof, in order to record what the RX antenna could “see.”
In addition, the RX and TX trajectories were recorded with accuracy at dm level, corresponding to the envisaged ns delay resolution. If for instance the LOS signal path is blocked, both positions are required to determine the delay of the hypothetical LOS path. For the van this was accomplished by storing data from wheel sensors, GPS and an inertial measurement unit. For the Zeppelin the situation was more difficult. Because of the high transmit power and being close to the GPS L1 center frequency the acquisition and tracking of GPS signals were impossible in a range of about m around the TX antenna. Hence, for a coarse positioning, an image of the airship was transmitted from the camera positioning station on the ground to the Zeppelin’s pilots, allowing them to hold the horizontal position within a range of about 30 meters. The vertical position could be kept more precise by use of the on board variometer and altimeter. The remaining movements were in the order of m/s and were logged on the so called Doppler station by means of recording the frequency changes of a GHz beacon transmitted from the Zeppelin. This allowed us to correct the influence of the Zeppelin movement on the channel measurement data in postprocessing.
2.2. Channel Sounding
For the measurement a commercial channel sounder of type RUSK [17] was used. This sounder transmits a measurement signal with period , generated from complex frequency domain samples via an inverse DFT (discrete Fourier transform) [18]: which filled the complete measurement bandwidth with equally spaced sine waves as depicted in Figure 2. Amplitudes and phases of the sine waves were adjusted to minimize the crest factor in order to have low distortion in the power amplifier of the transmitter [19].
To obtain the CIR from the received signal the receiver performs a DFT. Because of the periodic character of the signals and the synchronization between TX and RX, which was established by the atomic clocks, the momentary CIR can be estimated by cross-correlation of the received signal and the locally duplicated transmit signal over one period in the frequency domain where is the DFT of the received signal including additive noise [18]. is the estimated channel frequency response. Under the constraints of an additive white Gaussian noise signal, (2) represents the maximum-likelihood estimation of the Fourier transformed CIR [20].
The sounder settings for the measurement campaign are summarized in Table 1. Given the rectangular transmit signal spectrum, the time domain equivalent is a sinc function denoted as , obtained by Fourier transformation. For the bandwidth MHz the sinc function has zero crossings at with , thus the delay resolution or time between two uncorrelated signal samples in the measured CIR is
The repetition rate of CIR recordings, often denoted as snapshot rate , was chosen to avoid frequency ambiguities. The Nyquist theorem has to be fulfilled; that is, considering positive and negative Doppler the RX velocity during the measurements must not exceed
In the urban environment the RX velocity did not exceed km/h, allowing for extra margin, which was required as moving reflectors were present.
2.3. Antenna
The used RX antenna had a characteristic typical for satellite reception in land mobile applications, to be specific a hemispherical pattern and a right hand circular polarization (RHCP). A side view of the antenna is shown at the bottom of Figure 3. During the car measurement the antenna was mounted on the van’s roof rack (see Figure 1). A similar ground plate with the antenna in the center (see small image in the top left corner of Figure 3) was used on DLR’s rotational antenna test facility to measure the actual antenna pattern.
In Figure 3 the normalized copolar (RHCP) and cross polar (LHCP) reception for an azimuth are plotted at the measurement center frequency of GHz. On the -axis and correspond to the horizon and is the direction to the zenith. Being very typical for circular polarized antennas the cross polar suppression (see green curve) goes to zero for low elevation angles near the horizon. That means multipath suppression is ineffective at these low elevations. Ideally a reflected signal changes its polarization at the reflection point and so these signals would be nicely suppressed if they arrive from high elevations. In reality the reflectors are at very low elevations seen from the receive antenna. That is why cross polar suppression very often fails.
For the pedestrian measurements the antenna was carried by a person walking on the sideway as shown in Figure 4.
2.4. Scenarios
The urban measurements were conducted in the city of Munich with its population of 1.3 million, whereas the suburban measurements took place in Fürstenfeldbruck (34 thousand inhabitants). An important goal for the new channel model was to provide full satellite elevation and azimuth dependency. To reach this goal, measurement tracks were selected for the different channel types. Then for each type repeated drives and walks, respectively, along these track were measured with the Zeppelin at different altitudes and distances to realize different elevation angles between and . In total, 40 runs each lasting for 10 to 20 minutes were recorded.
3. Super Resolution Processing
According to the measurement bandwidth the delay resolution in the recorded channel impulse responses is ns. To estimate the precise delays of echoes we applied a super resolution algorithm.
For the estimation of the delays of several superimposed replicas of a known signal the maximum likelihood (ML) principle provides an efficient estimation but requires solving a multidimensional optimization problem [21]. One method to avoid this optimization is to calculate a suboptimal estimation like in MUSIC (multiple signal classification) based algorithms, for example, in [22]. Another strategy is to exploit the features of the considered signal to reduce complexity. Both approaches are combined in ESPRIT (estimation of signal parameters via rotational invariance techniques). ESPRIT is a frequency estimation technique which is based on a harmonic model; that means that the considered signal consists of complex exponentials in noise. Its ability to resolve complex exponentials closely spaced in frequency has led to the name super resolution [23].
In practice the performance degrades the closer the exponentials are spaced. For our analysis we applied the algorithm described in [24], which uses a simple interleaving technique that significantly improves the frequency estimation accuracy. Analogous to frequency estimation this algorithm can also be used to estimate specular components in the CIR from the measured channel transfer function. The achievable resolution is only limited by the signal to noise ratio and a remaining measurement device calibration error.
Super resolution processing in this application can be considered to be extrapolation of the measured channel power spectrum density; that is, the delay resolution of the channel impulse response is increased. Under the modeling assumption of distinct echo paths the CIR can be written as sum of delta functions. The Fourier transformation into the frequency domain gives
That means is periodic in this case. An echo with delay ns would transform into a sine wave with a period of exactly the measurement bandwidth of MHz in the frequency domain, as demonstrated in [25]. For higher delay resolution it is necessary to estimate contributions of sine waves in with lower frequencies, thus only sections of the sines periods are within the measurement bandwidth. Due to noise in the measurement data the frequency estimation of these sine waves becomes more difficult when the present echoes are spaced closer in delay.
As an example a single measured CIR is shown in Figure 6 together with the super resolution processing result. The red “+” are the measured samples and the blue lines with bubbles on top are the estimated distinct echoes estimated by the ESPRIT algorithm.
In Figure 7 a series of 500 consecutive recorded CIR snapshots is shown. For this plot the measured CIR samples with phases at delays were interpolated by with MHz. Being analogous, the super resolution data was interpolated with MHz according to an expected accuracy of ns in the ESPRIT result [24].
The power in Figure 7 is color-coded in dB and normalized to the unobstructed LOS signal. Like focusing with a zooming lens, the elliptical area shows the enhanced sharpness due to the super resolution processing. In contrast to the raw measurement data clear traces become visible. Each trace corresponds to a single echo and shows how its delay changed over time. In the beginning of this measured sequence the van was at a stop and thus the relative delays were quite constant. After about 150 snapshots the van started to accelerate and the channel began to change rapidly as many short existing echoes with changing relative delays and changing powers indicate.
4. Echo Signal Analysis
The super resolution processing result of the high resolution measurement data allows a deep “view” into the multipath propagation channel. The found traces verify the assumption of distinct echo paths. Moreover, the following analysis will show that most echo signals originate from small structures and objects in the near environment of the receiver and that the delay and Doppler frequency evolutions are mainly caused by the moving receiver.
4.1. Echo Detection
We are now going to apply a bounded multidimensional gradient search algorithm in order to automatically detect the start and stop times of echoes and to track their delay evolutions.
Starting from the sample with maximum power in the ESPRIT result, each echo’s delay is tracked in forward and backward direction of time . From a current sample at time all samples which are less than in the future or past, respectively, are analyzed. For all these samples with indices the weighted sum of the following four parameters is calculated as follows:(i)the delay distance with its weight ,(ii)the power distance with its weight ,(iii)the time distance with its weight ,(iv)the delay rate distance with weight ,
where , , , and were determined empirically. Note that, because the delay changes of echoes are increasing with higher speeds, the delay weight was implemented depending on the momentary receiver speed .
In the following step the sample with minimum is identified as the next point in the trace of the echo. The forward and backward detection end, if an upper bound for the multidimensional distance is exceeded: min, with . Due to the upper bound on the weighted sum as a stop criteria, the search range on the time axis for echo samples which are analyzed in a single detection step is given by .
The result of this echo detection (black lines in Figure 7) provides information and statistics for echo life cycles, delay evolutions, and the number of coexisting echoes.
4.2. Complex Echo Signals
In the next step the complex echo signals are extracted from the measurement data, exploiting the result of the echo detection algorithm with the delay evolution of each single echo. From the measured channel impulse responses with sample amplitudes and phases at sampled delays , the complex echo signal samples can be calculated by interpolation. For the rectangular band limited measurement signal spectrum with MHz the sinc function results as the interpolation function according to the sampling theorem [26], which allows the unique reconstruction of the continuous signal from samples of an ideal band limited signal where the Nyquist criteria is fulfilled:
4.3. Doppler Shift and Echo Bandwidth
Let be the complex signal samples at CIR snapshot times of a single echo determined in the previous subsection; then the signal spectrum can be calculated via a DFT where , is the number of echo signal samples, is the time of the first echo signal sample, is the CIR snapshot rate, and
In mobile scenarios the spectrum of each echo signal may significantly change over time. Thus we calculated so called spectrum snapshots at time instances from the echo signal samples in the interval . We chose as a trade of between frequency resolution and averaging effect.
Figure 8 shows a spectrum snapshots series calculated from a single echo signal. In this plot the power spectral magnitude is color-coded in dB. For this example echo, it can be observed that the Doppler shift falls from Hz to Hz. The reason is that at this moment of the measurements the receiver was just approaching and finally passing by a reflector that caused this echo. Note that the actual bandwidth of the echo signal is only a few Hz without significant changes over the echo’s life cycle of about s.
In the next step the Doppler shift evolution over time of an echo can be determined by searching for the maxima in the interpolated spectrum snapshots. The result is shown as black line for the echo example in Figure 8.
While the time variant Doppler shift was found to be introduced by the receiver movement, the additional phase variations in the echoes were caused, for example, when a reflection point moved over a rough surface or by the movement of branches and leaves. But also effects introduced by the “receiver” caused a certain spectral bandwidth, for example, vibrations of the van’s chassis on bumpy roads and swaying in curves and when breaking, or the regular antenna movements of a walking pedestrian could be observed in our measurement data and were verified by the on board video recordings. Analyzing all echo signals from various scenarios, the distribution of echo bandwidths (shown in Figure 9) was found to be uncorrelated to the receiver velocity and independent of the satellite elevation.
4.4. Fading
In narrow band channel models the Rice factor defines the ratio between LOS (constant) power and the (fading) power of all echo signals [26]. Similarly, we observed fading on each single echo signal extracted from our measurement data, resulting from the interaction of all electromagnetic waves reflected from one object. The observed fading characteristic can be described the same way in this subsection of the channel with where and are the mean and the standard deviation, respectively, of the echo signal’s envelope .
Analogous to the echo bandwidth, the Rice factor distribution shown in Figure 9 was found to be independent of the satellite elevation.
4.5. Reflector Location
In the single antenna case of the here presented measurements, the AOA of echo signals is a so far missing information. However, the extracted delay and Doppler trends of echoes, as well as the recorded receiver movement, can be exploited now to estimate the points of reflection in space. The basic assumption for the following analysis is first order reflections, which means an echo path passes only one object.
Assuming a stationary reflector, the relation between the Doppler shift of the received echo signal and the relative speed of the receiver to the reflector is defined by
The relative speed of the receiver results from the projection of the receiver speed vector onto the unity vector from the receiver to the reflector : with the angle between and . From (13) and (12) the angle can be calculated for a given receiver speed and a given Doppler of the received echo by
The situation is displayed in Figure 10 with the -axis being the heading direction of the receiver on board of the van. Seen from the moving receiver, the reflection point R has to lie on the surface of a cone according to (13), which is a first condition to be fulfilled. This cone has an opening angle of and opens in the driving direction for a positive Doppler or behind the car for a negative Doppler.
Let be a coordinate system with origin O at the receiver antenna location, with being the heading direction of the van, the lateral axis positive to the left, and the vertical axis; the unity vector to the reflection point on the cones surface is defined by where is the reflectors angular position on the cone (left hand positive, positive -axis ). In the same coordinate system the direction to the known transmit antenna position Z on the airship is with the azimuth direction (left hand positive, front ) and the elevation to the airship (zenith ) as shown in detail in Figure 10.
The second condition for the reflector position is related to the excess delay of the echo signal. The known distance between the two antennas plus has to be equal to the echo path length: where is the distance between RX and reflector. Geometrically this condition results in an ellipsoid with the receive and transmit antenna in its focus points. The sum of the distances from any point on the ellipsoid to the foci is the echo path length, which is equal to the major axis of the ellipsoid, and the linear eccentricity [27] is given by .
By inserting (14), (15), and (16), the two remaining unknowns in (17) are and the angle which define the reflector position. The solution is the intersection line of the cone with the ellipsoid (red dash-dotted line in Figure 10).
As a third condition we now look at the reasonable -coordinate range for the reflection point R on this intersection. In the city center of Munich the average building height is about m. Taking the height of the cars antenna platform of m into account, the range is . According to (15) the -coordinate of the reflection point is given by . For a hypothesis of the distance is therefore and together with (17) the solution for the reflection point can be calculated numerically. Since facades of houses and most other structures in urban environment are vertically aligned, we further assume that the elevation angle of the echo path does not change at the reflection point; that is, sin shall be fulfilled.
We found three possible cases.(i)One solution: for most of the detected echoes a distinct reflector location can be determined.(ii)Two solutions: when the satellite incidence azimuth is almost parallel to the van’s heading, the reflection point may be on the left or right side of the street. We chose the one with lower , which lies on the more illuminated side, to be more likely.(iii)No solution: the condition sin can only be fulfilled for . This means the reflection comes from a nonvertical structure. Therefore we set modeling a reflection on the roof of a building.
A snapshot from the car measurement in the center of Munich in Figure 11 illustrates the result of the above analysis. In this aerial photograph the receiver track (green line) with the power of the direct path (size of the green bubble) and several echoes (red bubbles) are displayed. In this situation the direct path signal was about dB attenuated by a tree on the left side of the van. At the same time two strong echo signals were received from the house just to the right, and several more echoes with lower powers originated from nearby houses, trees, and other objects in the street.
Satellite Elevation Dependency. We define polar coordinates , where is the norm of the projection of onto the horizontal plane and is the incidence azimuth of an echo signal relative to the heading direction of the receiver (-axis). For each of the measurement runs (with nearly constant elevation) we determined that is, the measured reflector location likelihood averaged over all satellite azimuths weighted by the azimuth probability on the specific measurement track. An example of this distribution is shown for the urban city car environment at elevation in Figure 12. It is notable that the likelihood of reflectors in front or in the back of the car (along the -axis) is very low, while the majority of reflections is received from both sides of the street. The analysis of data at other elevation angles reveals that the azimuthal area in the front and back of the car with low reflector probability increases almost proportionally to the elevation angle. Thus, the reflecting area is cut back to two narrow zones at the sides of the street, when the Zeppelin is at the zenith. On the other hand reflections from a wide zone along both sides of the street were received, when the elevation angle was below .
4.6. Relative Reflector Azimuth
We define the angle as the azimuthal angle of arrival of echoes relative to the satellite azimuth. The discrete probability density function can be determined from the processed measurement data. Similar to the horizontal reflector position distribution we are taking into account that the satellite azimuth was not equally distributed during the measurement; therefore the likelihoods are weighted by the occurrence probability of the present satellite azimuth:
As an example this distribution is plotted in Figure 13 for the car measurements in the urban city center environment at elevation. The plot shows that the likelihood maximizes around . This maximum can also be observed for the other measured elevations. Consequently for a car in an urban environment a satellite on the left side most likely causes reflections on the right side of the street. More surprisingly the data show that even when the satellite azimuth is to the front left position of the car, most of the reflections come from the right side behind the car and not from the front right, where the specular reflection point would be on a flat house wall.
The reason for this very interesting result is double reflections on corners illustrated in Figure 14, which occur very often in cities at, for example, projections or ledges on walls. The main characteristic of such a reflection is that a wave coming from the satellite is reflected back into the satellites direction in the horizontal plane. This result is comprehensible, because in contrast to a single reflected signal which becomes cross polarized for circular polarized transmit signal, a double reflected signal is turned back to copolarization; thus it will not experience the typical cross polar suppression at the receiver antenna, and its magnitude at the receiver front end can be even higher than from single reflected signals. We want to point out here that, due to the lack of scenery details in ray tracing simulation based channel models, such models are not able to reproduce this important effect, which has a strong impact on the incidence angle distribution of echoes.
As stated in the beginning of this section the underlying assumption for the reflector position estimation is first order echo paths which pass only one object. For instance in case of a double reflected path, a fictitious single reflection point is determined which models the delay and Doppler trend observed in the measurement. The so-gained statistical single reflector position distribution can still be used in the model to reproduce the real channel conditions with respect to the Doppler and delay trends of echo signals. Moreover, for the double reflections on corners which we identified as major multipath source together with the single bounce, the resulting error in estimating the incidence angle is negligibly small.
5. Modelling
We propose a novel Wideband MIMO Satellite (WMS) channel model that reproduces multipath propagation effects in accurate detail by taking correlations among echoes from single reflectors into account. The description starts with the MISO (multiple input single output) case, where we consider a single receive antenna, and will be expanded at the end to MIMO configuration.
5.1. Model Concept
The basic model concept is a combination of deterministic and statistical modeling. Several channel characteristics were identified to be dependent on the receiver position and movement, respectively. Thus, for instance the Doppler offset and delay trends of echo signals or the shadowing effect of houses on the direct signal can be calculated with low computational effort. On the other hand, the changing number of coexisting echoes, their bandwidths, amplitude variations, and life cycles are most realistically modeled by the statistics determined from the measurements.
The block diagram in Figure 15 gives an overview of the new model. The time variant speed and heading of the receiver, and the satellite elevations and azimuths are required input to a motion model that drives an artificial scenery, where the receiver moves along a straight street. The origin of the coordinate system is at the RX starting position and the -axis points into the RX heading direction along the street. In the motion model the RX position on the -axis is calculated from the user speed by integration and the satellite azimuths are transformed into the RX heading direction by that is, are the satellite azimuths relative to the user heading. With this approach a turn in the road can be modeled by simply changing the user heading, which leads to changes in .
Beside the street, the artificial scenery consists of house fronts, trees, and lamp posts affecting the direct signal paths from the satellites. In order to model different conditions in terms of visibility, the size and position of these elements are parameterizable.
For the multipath components a time variant and RX position dependent number of virtual reflectors are created for each of the satellites in the artificial scenery (see the illustration example in Figure 23). The reflector positions in the scenery are chosen according to the likelihood distributions determined from the measurements. The powers, bandwidths, Rice factors, and life spans of the resulting echoes are taken from the measured statistics too. As the receiver moves through the artificial scenery, the amplitudes (path index and satellite index ), phases , and relative delays of the echoes change according to the position and time dependent processes that we identified in Section 4, which we modeled both, deterministically and statistically. The model output is a series of complex time-variant channel impulse responses for each satellite link.
5.2. Direct Path Modeling
We define the direct path as a representation of the LOS transmission in the unobstructed case or as the attenuated and delayed path (with respect to LOS) due to diffraction in case of LOS blockage. In the measurements three major types of obstacles were identified that influence the direct signal reception: house fronts, lamp posts, and trees.
5.2.1. House Fronts
We observed characteristic power decays when the receiver approached obstacles, entering their shadowed area. In Figure 16 the power of the direct signal is plotted for a short part of the car measurements in the urban city center environment with the Zeppelin at elevation. In this situation the van drove towards a crossroad with a continuous row of houses ahead (house height m) and the Zeppelin was behind these houses. When approaching the crossroad the distance to the house front decreased and the van entered the shadowed area (right of the dotted line). From this instance on the signal is strongly attenuated.
We model this behavior with the so-called “knife edge model” [28]. In this model it is assumed that a planar wave hits a half-sided infinitively large conducting plate. The focal length is given by where and are the horizontal distances from the house front with height to the receiver and to the satellite, respectively. With the clearance distance the normalized Fresnel variable (see [28]) can be calculated by where is the wavelength. Positive values of correspond to positive values of ; that is, the receiver is in the illuminated region, whereas a negative corresponds to the geometrically shadowed region. In case of knife edge diffraction at a half-sided infinitively large conducting plate the complex diffraction coefficient at the normalized position is given by [28] where is the result of the Fresnel integral
The magnitude square hereby represents the intensity of the diffracted field relative to the intensity of the unobstructed field
The calculated signal power variation from this model is displayed in the left part of Figure 16 too. The comparison to the measurement proves the appropriateness.
Besides the diffraction on roofs, we also consider diffraction along the house edges. This is necessary to correctly model the appearing and disappearing direct signal when the receiver is passing for instance a crossroad or a gap in a row of houses.
From the known geometry in the artificial scenery the relative delays of diffracted signals can be calculated. In the example model output in Figure 23 the LOS path of satellite 1 (red) is blocked very close to the roof edge of a house front in the artificial scenery. The corresponding CIR in Figure 24 shows two diffracted signals with very short relative delay (marked as direct paths). One is diffracted at the roof (dB) and the other one at the closer house front edge (dB).
5.2.2. Lamp Posts and Tree Trunks
Lamp posts have a significant effect on the direct path too. On the right side in Figure 16 the power of the direct signal is plotted for a situation where the van passed a lamp post with a diameter of cm. The -coordinate denotes again the RX position along the street, with at the point where the LOS passes through the lamp posts center. The van’s horizontal distance was approximately m at this point. The measurement shows that the power of the direct signal begins to oscillate, goes down in the direct shadow of the lamp post, and increases again oscillating as the van moves away.
We model this effect with a “double knife edge model” where two overlapping plates are alternately present. One spans from and the other one from . Both edges are calculated separately and then added coherently. Thus, the complex diffraction coefficient is and the resulting intensity of the diffracted field is
Figure 16 shows the very good match of model and measurement. Note also the mismatch for as we approached the next obstacle in the measurement.
The same effect was observed at low elevation angles when passing tall trees, where the tree trunk caused the power variation. Therefore tree trunks can be modeled in the same way by this receiver position dependent process. In contrast to house fronts, the shadowed region of lamp posts and tree trunks is very small. Consequently, the relative delays of the diffracted signals are so small that we neglect them in the WMS channel model.
5.2.3. Tree Top Model
Analysis of our measurement data showed that the signal attenuation caused by trees mainly depends on the length of the direct path through a tree top. In addition we observed power variations, which were caused by the inhomogeneity of the tree top. This is because in L-band the wavelength is in the order of the size of branches and leaves and the power variations occur as the receiver moves through the nonuniform shadow of the tree.
In both, [29, 30] the attenuation through trees is determined from a complex tree scattering model, where branches are modeled as dielectric cylinders and leaves are modeled as discs with complex permittivity. These deterministic models are accurate but they require large computational resources.
Following our general hybrid approach for the WMS channel model, we model a tree top as one dielectric cylinder and we add a stochastic process to account for the variations in the shadow of the tree. The cylindric shape was selected to fit the trees along the measurement track (see Figure 5). Parameterized by the attenuation constant , the diameter, and the height of the cylinder, different kinds of trees can be placed in the artificial scenery of the WMS channel model to create specific scenarios.
(a)
(b)
The attenuation through the cylinder , which depends on the receivers position , can be calculated by where the direct path length within the tree top is geometrically calculated given the constellation in the artificial scenery.
The stochastic process is derived from analysis of the fluctuations in the measured direct path signals (caused by the inhomogeneity in the tree top). The spatial spectrum was found to have a Gaussian shape, which can be parameterized by the double sided dB bandwidth with unit 1/m. Furthermore the signal envelope was found to be Rician distributed, which we parametrize by the Rice factor : where and are the mean and the standard deviation, respectively, of the envelope.
We calculate the stochastic process by where represents the constant part of the Rician process with a randomly chosen phase uniformly distributed in . The stochastic process in (32), given by the sum of the constant and the fading part , is normalized by the factor . The fading part is derived from the normalized Gaussian power spectrum where the standard deviation is given by
As it is proved in [31], the fading part of can be represented by a sum of carriers with Gaussian distributed frequencies (according to the observed spectrum ) and uniformly distributed initial phases :
In the WMS channel model implementation we use a very efficient realization of (35) using “controlled randomness” [32] for improved performance versus complexity.
Finally the total field attenuation coefficient through the cylindric tree top model, where the Rician fading process overlays the deterministic attenuation of the cylinder, is calculated by
5.3. Echo Path Modeling
So far the correlation between direct paths from different satellites (or different antennas on a hovering platform) is well accounted for by the common artificial scenery in our WMS model. An open question that we want to answer in this section is how correlated the echo paths on different links are, and how can that be modeled efficiently?
5.3.1. Echo Life Distance
In the measurement data the channel appears rapidly changing. Many echoes disappear and new ones appear (see Figure 7). This process was found to be highly correlated to the receiver speed. When the car stopped most echoes remained with slow fluctuations. Relevant to our motion dependent model is the life distance of an echo, which we define as the distance the receiver is traveling from appearance of an echo until it disappears.
In Section 4.1 the life cycles of echoes were determined. The life distance is calculated by where is the time when the echo was first observed and is the recorded speed of the van. Figure 17 shows the cumulative distribution function of the life distances obtained from the measurement in the urban city center environment. It can be seen that most echoes exist along a motion path below m. They originate from objects like, for example, window frames or traffic signs, which are “visible” only for a specific narrow angle as illustrated in Figure 18.
The opposite way around when looking at the illumination side, we can conclude that echoes from this reflector can only be present simultaneously on the MISO channel for satellites whose incidence waves are separated by less than . Therefore we call the correlation angle which is specific for each reflector and depends on and on the distance between RX and the reflector . In the WMS model this correlation is implemented in the following way: whenever a new reflector is created on one of the multiple input channels, copies will be initialized on those channels which are separated less than , but their presence along the RX trajectory is shifted the more they are uncorrelated. Because of the pure horizontal movement of the RX, we consider only the azimuthal separation to be relevant in this modeling approach. For instance if two satellites are separated by in azimuth, then the shift along the -axis is . As a result echoes from this reflector are present simultaneously on both channels along of the receiver trajectory.
5.3.2. Number of Coexisting Echoes
In case of GNSS, the satellite constellations are designed to provide the best service by a steady high degree of separation of satellites in the sky [33]. Thus most single echoes and their corresponding reflecting objects will be independent in the MISO case. Nevertheless one can observe that clusters of reflections originating from a house front do look similar, even if the incidence angle changes for several . An important parameter which allows us to model this observation is the number of coexisting echoes. It is so important because the performance of multipath estimation algorithms strongly depends on the assumed number of paths [34]. As an example the changing number of coexisting echoes determined in Section 4.1 is plotted for a full measurement run in Figure 19.
We note a relatively high number of echoes (up to ) at the same time. Besides the mean number of coexisting echoes we incorporate the process of increasing and decreasing number of echoes in the WMS channel model. A Fourier analysis reveals two processes: a narrow band process with high power and a lower powered wide band process . The narrow band process reflects the changing multipath conditions in different streets with different house types, while the high frequent variations are caused by changes in facades and by passing smaller objects like trees or traffic signs.
In the WMS channel model a number of echo generator (NEG) controls the number of echoes accordingly for each satellite channel. In order to model the mentioned correlation between different channels we use the same realization of the narrow band process on all channels scaled by a factor which depends on the actual elevation, while the mean number of echoes and the wide band process are synthesized individually for each link.
Analogous to the other WMS model parameters, the number of coexisting echoes is modeled elevation dependent on each environment. Therefore the parameters for the wide band process, as well as the mean number of coexisting echoes and the scaling factor , were determined for all measured elevations between and for each environment. The exact parameter values are documented in Recommendation ITU-R P.681-7 [35] and Report ITU-R P.2145 [36].
5.4. Statistical Echo Signal Generation
While the NEG controls the partly correlated number of echoes on the multiple satellite channels, in general all the parameters describing a single echo signal are chosen independently of statistics as indicated in the block diagram in Figure 15—except for reflector copies (same position, echo power, and initial phase) in case of closely spaced satellites, with life cycles shifted along . Despite this fact, the representation of statistical reflectors as source of echo signals in our artificial scenery will allow taking their correlated deterministic phase and delay changes into account, as we will explain in this section.
5.4.1. Reflector Location
In (19) and (20), respectively, the horizontal reflector position distribution and the relative reflector azimuth distribution were determined for the measured elevation angles and , . Both distributions are independent of the actual satellite azimuth (relative to the van’s heading) as they were calculated by averaging over all possible azimuths.
In the absence of sufficient amount of data to prove or disprove their statistical independence, we approximate in the next step the common distribution for a given satellite elevation and azimuth by the product given the definition of the relative reflector angle .
The result of this operation is shown in Figure 20 for and . Before the multiplication both distributions were calculated for from the measured distributions at and elevation using linear interpolation. In the plot the thought receiver position is at the origin and the van’s heading direction is along the -coordinate (red arrow). The satellite is in front of the van slightly to the left and the most likely reflector positions are on the right side of the street just behind the van. For comparison see the characteristics of the unconditional distribution in Figure 12.
In the WMS channel model new reflectors are created at horizontal positions according to this distribution depending on the satellites elevation and azimuth. For the selection process our implementation uses the rejection method [37]. The complete set of used statistics and parameters are documented in [35, 36].
Next the dedicated vertical position of each reflector is calculated applying the same rules as in the analysis step in Section 4.5: buildings with vertical walls dominate in urban and suburban environments. In case of a specular reflection on a wall, but also for the corner reflection, this means that the elevation angle of a ray does not change at the reflection point and therefore equals the satellite elevation . A virtual house front with height can be thought at the horizontal reflector position at the distance from the van. Thus, can be calculated geometrically by
The house height is chosen from the same statistics used to generate the artificial scenery. If would exceed of the virtual building, it is reasonable to assume a scatterer at the roof rather than a specular reflection. In this case is limited to , where is the receive antenna height above ground. Therefore the following rule is applied to calculate : which finally leads to the 3D reflector position.
5.4.2. Mean Power of Echo Signals
For each elevation and environment we determined the average power distribution and the power variance distribution of the measured echo signals. In general terms these statistics show (see example of in Figure 21) that the most powerful echoes originate from reflectors on both sides of the street close to the receiver. With increasing distance of the reflector, the average power decreases. But especially for elevations below a significant contribution from objects up to 100 meters in front or in the back of the receiver close to the street was observed.
The individual mean power of each generated echo signal in the WMS model is drawn from a Gaussian distribution parameterized by ( and ) according to its horizontal position and the momentary satellite elevation . Again the elevation dependency is realized by linear interpolation from the measured statistics. All modeled distributions are documented in [35, 36].
5.4.3. Time Series Characteristic of Echo Signals
When a car drives through a multipath environment, the receiver is moving through a quasistationary field radiated by reflectors. We are going to model this process in a deterministic way depending on the receiver position. But we also observed during the measurements that fading on the channel does not completely stop even when the car comes to a stop. This is because the channel is changing, for example, due to trees in the wind, other cars, and pedestrians. Furthermore, in our analysis we did not find a correlation between this fading process and the receiver speed. Thus, we will model this stochastic process time dependent only.
In the WMS channel model the echo signals are generated combining these deterministic and stochastic processes as illustrated in Figure 22. For the deterministic part the geometrical reflector representation is used. A reflector is initialized at a random position according to the distribution given by (38). When the receiver moves along the street in the artificial scenery, the path excess delay and the Doppler phase of the echo are calculated geometrically. The excess delay is given by where and are the vectors from the receiver to the reflection point and to the satellite, respectively. The scalar product divided by is the norm of the projection of in . The phase change in the reflected signal only depends on and the wavelength and is calculated by
This causes the deterministic process of the phase variation to be motion dependent only.
In order to model the stochastic process, we analyzed the spectrum and amplitude distribution of the echoes. The spectra were found to have a Gaussian shape with a typical double sided dB bandwidth in the range of a few Hertz. In Figure 22 a filtered noise signal represents this stochastic process. The observed Rician amplitude variation of a single echo is modeled by the Rice factor , which is defined as the ratio between the constant power of the echo and the power of its fading process . For each newly initialized reflector a bandwidth and a Rice factor are drawn from the distributions shown in Figure 9.
The stochastic part of the echo signal with mean power can be calculated by where represents the constant signal of the Rician process with a randomly chosen phase , uniformly distributed in the interval . The fading signal is derived from the normalized Gaussian power spectrum where the standard deviation is given by
The process in (43) given by the sum of the constant and the fading signal is normalized by the factor . In the WMS channel model the fading signal is realized as a sum of carriers (similar to (35)) with Gaussian distributed frequencies and uniformly distributed initial phases :
The complete echo signal is than calculated by where the deterministic Doppler phase due to the receiver movement and a random initial echo phase are added to the stochastic signal.
5.5. Multiple Receive Antennas
So far we have considered multiple TX but only a single RX antenna. Since our model approach is based on a movement along the -coordinate, we can simply recalculate the generated CIRs for different positions at a given time instance . In this way, for example, the diversity gain for dual antenna receivers on a car (one antenna in the front and one at the rear) can be evaluated. But we can also apply different movement models to each RX antenna to investigate scenarios where cars in a convoy change their relative distances. Such a scenario is particularly interesting to investigate the performance of new safety systems (e.g., using direct car-to-car communication) where precise relative position information is required.
On the other hand we are not free in placing RX antennas on different lanes on the street. Our measurement data does not provide the channel cross-correlation information for such a case, because we measured with only one antenna along one track. Nevertheless, Figure 17 indicates that a full correlation can be assumed in the first approximation if the antenna spacing is below cm, which is about half the wave length at our measurement frequency. Based on this approximation we can generate MIMO channels for a 2-by-2 GNSS antenna array like the one presented in [38] with full multipath correlation in - and -direction and the measured and modeled decorrelation along the -direction. Of course when calculating echo and diffracted direct signals, the delays and phases change depending on the exact position of an RX antenna element according to the geometry in the artificial scenery of the WMS model. That is, phase shifts in the generated MIMO model output reflect the incidence wave angles of the signal paths.
6. Generation of CIR Time Series
We implemented the WMS channel model in MATLAB. It can be downloaded from [39]. Configuration files are provided which allow specifying the user type (car or pedestrian) and the environment (urban or suburban) including parameters for the generation of objects in the artificial scenery. Movement models for the RX antennas can be applied and the incidence angles of the (unlimited) number of satellite waves can be varied throughout a run.
For reasons of clarity we show in the following example output a MISO configuration with a constellation of 4 satellites. The CIRs were generated at a rate of Hz on an 8-core CPU in real time. A snapshot of the channel state in an urban environment with a receiver in a car is illustrated within the artificial scenery in Figure 23. Four different LOS paths are visualized together with their echoes—respectively, the reflector locations where they originate from and their momentary powers (size of the spheres)—in the corresponding colors: channel number 1: magenta, channel number 2: red, channel number 3: turquoise, and channel number 4: blue.
The classic view on the CIRs for the same snapshot can be seen in Figure 24. On each channel the direct rays are plotted in red with a big circle marker and echo paths are blue with small markers. In this situation there are two diffracted paths on channel number 1 from the roof and the side edge of a house front as the LOS is blocked. Similarly the first path on channel number 4 is attenuated and slightly delayed with respect to LOS due to shadowing by a house front and in addition by a tree, which is typical for this low elevation of . Because of the chosen building height distribution replicating the measurement area in the city center of Munich, the shadowing effect of houses is diminishing for elevations above . In fact we see in this example that on channel number 2 the attenuation of the direct path is the result of an occasional single tree in the middle of the street, whereas satellite number 3 at is unobstructed. According to our findings this snapshot also resamples the decreasing delay spread for higher elevations and it shows that the mean number of echo signals is the highest at elevations where the house facades and objects on the side of the streets are well illuminated—here at about . For lower elevations the houses on the opposite side of the streets cause shadow even for the potential reflectors, whereas for higher elevations the echoes from vertical reflectors can only be seen on the side walk but not any more in the middle of the street.
Both types of channel state illustration can be watched throughout a run of the implemented channel model. When the receiver is at a halt, the phases of the rays in Figure 24 rotate slowly and the sizes of the spheres in Figure 23 fluctuate moderately. As soon as the velocity rises, the motion dependent processes cause fast changes in both plots.
The practical relevance of the ability to simulate such scenarios with our detailed model approach is demonstrated in [40]: simulations are performed based on the former SISO implementation of the model. Results for different receiver architectures, as well as for different signal structures, are presented that allow the assessment and comparison of the individual ranging performances. Important to note is the fact that large ranging errors with low occurrence probabilities are realistically reproduced, which is an important issue in safety critical applications. Further analysis of the simulation results also allowed for a better understanding of the mechanism that can cause loss of lock of GNSS signals in situations with disadvantageous changes in the relation between the direct signal and echoes, respectively, their amplitudes, phases, and delays.
7. Model Verification
For verification of the model implementation, we simulated routes that matched the actual measurement drives; that is, we entered the receiver heading and speed, as well as the incidence azimuth and elevation direction recorded during the measurement, and compared the output statistics of these SISO simulations to the measured statistics.
For each simulation and the corresponding measurement run, we calculated , that is, the total number of echo signal samples normalized to the snapshot rate , which was Hz in the simulation via Hz during the measurement. The example in Figure 25 shows the comparison for different elevation angles in the suburban car environment. It verifies the implemented generation of the number of coexisting echoes and the modeling of the echo life distance distribution. The general slight increase in the number of echoes for increasing elevation is reproduced.
Next we look at the power delay profiles. In particular we compare the discrete probability density functions of the expected value of coexisting echoes from simulation and measurement, respectively. From these distributions in Figure 26 we can read out the occurrence probability of an echo within a delay bin of ns and a power bin of dB. Their surface integrals equal the corresponding mean number of coexisting echoes. As an example Figure 26 compares the suburban car environment for the lowest measured elevation of and a high elevation of . The plots show that the channel delay spread decreases with increasing elevation and confirm the statistical fit of model and measurement.
8. Conclusion
High measurement bandwidth and precise synchronization allowed us to gain a unique insight into the satellite propagation channel and to isolate single reflections from small structures. The spatial reflector distribution that we found, respectively their back scatter properties, clearly defer from assumptions of previous models and show limitations of physical simulators in complex environments. We developed a comprehensive model that takes several important correlations into account. Specifically the delay and phase trends of coexisting echoes and the path delays relative to the hypothetical line of sight in case of diffraction are essential to correctly model the effects that cause large positioning errors in GNSS applications in multipath environments. The correlation in the MIMO channels is modeled individually for each reflector depending on its spatial characteristics. A combination of deterministic and statistical modeling in time and space domain allows fast and realistic generation of CIR time series for specific user scenarios, for example, to investigate movements like stop and go. The model is targeted to support the development of new communication and navigation receivers that can mitigate performance degradation due to multipath and shadowing. The applicability ranges from GNSS positioning, pedestrian navigation, and collision avoidance in intelligent transport systems to communication links to hovering platforms, where high availability of service and high accuracy are demanded.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to thank Fernando Pérez Díaz and Thomas Jost for implementing the MISO extension of the model.