On the Feasibility of High Speed Railway mmWave Channels in Tunnel Scenario

1State Key Laboratory of Rail Traffic Control and Safety, Beijing Jiaotong University, Beijing 100044, China 2Beijing Engineering Research Center of High Speed Railway Broadband Mobile Communications, Beijing 100044, China 3Mobile Application Research Department, Electronics and Telecommunications Research Institute (ETRI), Daejeon 34129, Republic of Korea 4School of Electrical Engineering, Korea Advanced Institute of Science and Technology (KAIST), Daejeon 34141, Republic of Korea


Introduction
Mobile communication is one of the most essential and successful technology progress types, which becomes an indispensable part of more than 5 billion people [1].With the mobile data demand exponentially growing [2][3][4][5], wireless communication systems working at millimeter wave (mmWave) have attracted much more attention in many mobile communication scenarios.For mmWave, it potentially contains a large amount of spectrum resources for achieving multi-Giga bps data rate for wireless communication systems [6][7][8].
As a typical application scene of mobile communication, high speed railway (HSR) with speed of more than 300 km per hour challenges the constantly improved mobile communication systems [9].The most influential challenges for HSR mmWave communication system include frequent handovers, difficult signal processing as very high speed, and high penetration loss of signal from base-station to intrawagon user [10].However, for achieving multi-Giga bps data rate in high speed railway (HSR), the existing HSR dedicated communication systems are no longer satisfactory.For example, the maximum data rate for the Global System for Mobile Radio Communications for Railway (GSM-R) is less than 200 kbps.Even for the Long Term Evolution for railway (LTE-R), it cannot provide more than 100 Mbps data rate.Therefore, an interest group of High Rate Rail Communication (HRRC) has been established in the IEEE 802.15 working group for developing advanced mobile communication technologies, and it targets final achievement of very high data rate for HSR wireless communication system.Also, a distributed antenna system (DAS) [11] has been designed and it becomes a very promising communication system for HSR and metro system [10].However, all these desiderate full understanding of the channel and an adequate and reliable channel model [12].

Wireless Communications and Mobile Computing
Since channel sounding measurements are treated as a priority when studying channels, in the era of 3G and 4G communication system, the channel parameters and models are extracted mostly from massive measurement campaigns [13].But the measurement at mmWave band greatly challenges the existent sounding systems, especially for highdynamic channel measurement.The HSR mmWave channel measurement data is even rare [3].So, we propose an approach that enhances the capability of deterministic channel model for obtaining mmWave HSR channel.Then, some important channel characteristics are studied and modeled for guiding future HSR mmWave communication system design.
In this paper, we simulate an important wireless channel by an enhanced ray tracing simulator where the electromagnetic (EM) and scattering parameters of dominant materials in the scenario are measured and fitted with suitable models.The simulation configuration practically follows the communication system demands in [15].Meanwhile, a timeinterpolation method is employed for promoting the simulation capability for high-dynamic channels.Based on extensive simulations, the channel characteristics in an arched tunnel are completely exposed in two defined regions.Moreover, the channel feasibilities of using different antennas are investigated for determining the most suitable one.
The remainder of this paper is organized as follows: Section 2 describes principles of the proposed ray tracing tool with material measurement and time-interpolation.In Section 3, simulation scenario and system setups will be presented with more details.Section 4 demonstrates the feasibility of communication system with different antenna setups.Section 4 also illustrates the necessity of partition of regions in radio channel analysis.Then, characteristics of radio channel are provided in Section 5. Finally, Section 6 gives the conclusion.

Ray-Optical Based mmWave Channel Modeling for High-Dynamic Scenario
The wireless channel is simulated by utilizing ray-optical based deterministic propagation model, that is, ray tracing (RT) model.As it differs from other traditional RT, the enhanced RT used in this study has been calibrated in various frequencies and scenarios [16][17][18][19].We make a special measurement for extracting electromagnetic (EM) and scattering parameters of the most influential materials in the HSR tunnel.The extracted parameters will be implemented into RT simulator.Moreover, for investigating the detailed channel characteristics in such high-dynamic scenario, a well-studied RT time-interpolation algorithm is employed to expose the channel small-scale characteristics.
The RT, in this study, is developed based on ray-optical algorithm [20] at Technische Universität Braunschweig [21].It is three-dimension (3D) channel simulator which is performed in 3D digital map.For an integrated channel modeling, several radio propagation mechanisms are taken into consideration as line-of-sight (LOS) ray, reflection rays, and scattering rays.The LOS ray is known as free-space wave propagation between transmitter (TX) and receiver (RX).Its power is calculated by free-space path loss model in RT.The LOS path, if it exists, will dominate the received power in the RX.Then, the reflection which is defined as the incident angle and reflected angle of EM wave from a surface is the same.Its power is calculated by solving the well-known Fresnel equations and the related reflection point in the 3D digital map is obtained by applying the image method [20].Still, in the case of reflection, th order multiple reflections will be practically considered for capturing dominant power contributions of total received power.For the diffuse scattering, the power contribution of diffuse scattering is evaluated by adopting the "effective roughness" (ER) scattering models [22]; the ER models generally include Lambertian Model, Directive Model, and Backscattering Lobe Model for predicting different types of scattering patterns [22].Owing to the special structure of the tunnels where no obstacles are situated, the transmission through the materials and diffraction from edges of obstacles are omitted in the simulation.
Based on ray-optical principle and several propagation mechanism models, the final output of the RT is a timevariant channel impulse response (CIR) ℎ(, ) which is the sum of powers of all the determined rays [23]: where   (),   (), and   () denote amplitude, delay, and additional phase shift of th ray (totally () rays that have been found by RT kernel).
Then   () =   ()⋅ (2  ()+  ()) is defined as a complex coefficient; formula (1) can be rewritten as Then, (2) attached with the effects of wave polarization and antenna gains;   () can be further expressed by where  RX/TX, and  RX/TX, indicate the angle of arrival (AOA) or angle of departure (AOD) of the th multipath component.Therefore, the vectors  →  ( RX, ,  RX, ) and  →  ( RX, ,  RX, ) denote the complex polarization vectors [24] of antennas at RX and TX, respectively. RX, and  TX, indicate the additional antenna gains of the th multipath component in RX and TX.(⋅)  denotes the Hermitian transpose.  () depicts the complex channel polarization matrix which indicates the polarization shifting of the ray.  () comprises propagation loss and the phase shift according to the delay   () of the th multipath component.More information can be found in [21].
For investigating mmWave channel, it is still the challenge that the mmWave channel sounding measurements are very costly and time-consuming even for static scenarios.For measuring high-dynamic mmWave channel, it is even unrealistic to measure a mmWave channel sample at a normal HSR speed (around 300 Km/h in China).So the RT is a promising tool for revealing the channel characteristics in HSR scenarios.It is also a trade-off between the limitation of measurement and demand for the channel data.To enhance the RT capability of accurately predicting channels in HSR tunnel scenario, we perform two approaches, as mentioned in the initial part of this chapter, which are measuring the EM and scattering parameters of the most influential material and utilizing RT time-interpolation algorithm to extract the small-scale channel parameters.

EM and Scattering Parameters Acquisition of the Material.
The accurate energy calculations of the reflected and scattered rays in RT simulation depend on accurate EM (relative permittivity   in our study) and scattering parameters of materials as well as accurate 3D digital map of the scenario.Generally, RT simulation starts with determining the materials which constitute the scenario; then   and scattering parameters of these selected materials are acquired from checking literatures (e.g., ITU-R recommendations [25]), performing dedicated EM measurement [26], or derivation from channel measured data [27].However, apart from dedicated EM measurement, these methods are less than ideal as   and scattering parameters of a material vary in different scenarios and conditions (e.g., humidity of environment and density of material).Therefore, we measure the interferences of the material on the wave propagation at 32.5 GHz using a self-built testbed and estimate   and scattering parameters accordingly.And the material, in this study, is sulphoaluminate cement which is generally applied in HSR tunnels (cf. Figure 1(a)); unlike classic cerement utilized in buildings, it has advantages of quick-drying, high strength, compaction, and cementation which requires a dedicated measurement for studying its unknown   and scattering parameters.
As it is shown in Figure 1(b), the self-built testbed is on the basis of two high-accurate rotatable arms which could lead an accurate 2D scanning of the scattering of a material; the TX and RX are with horn antennas; and the Vector Network Analyzer (VNA) used in the measurement is manufactured by Keysight Corporation with model N5247A.The VNA measures  −  between RF-ports of two cables after the end-to-end calibration process.In the measurement, the material is hanged upside the rotating center by a rope for eliminating some unwanted interference, which is consistent with the function of the anechoic chamber.As it is illustrated in the right of Figure 1(a), the diffuse scattering data is obtained by rotating the RX while the TX is fixed; and when  Inci =  Scat , we get the reflection data.After measurement, the data are postprocessed with proper filter to attenuate some interferences.Then, we can estimate   of the material by the method of free-space measurement based partly on [28]; and the estimated   is verified and slightly tuned with the measured refection coefficient accordingly (cf. Figure 2(a)).Meanwhile, the scattering parameters can be estimated by the method similar to [22,26] but additionally with Simulated Annealing Algorithm [29] to obtain the better scattering parameters by automatic minimising of the gap between fitting results and measurement (cf. Figure 2(b)).It should be noted that the Directive Model as one of most important ER models was employed to fit the measured scattering radiation patterns; therefore frequency-dependent () and   () of Directive Model are the scattering parameters that should be extracted from estimation process [22,26].
As it is described above, a brick made of sulphoaluminate cement was measured; Figure 2 illustrates the estimation processes of relative permittivity (  ) and scattering parameters ( and   ) at 32.5 GHz. Figure 2(a) gives the comparison of reflection coefficients between fitting curve and measurement.  of the fitting curve is   = 3.47 − 0.15, where the mean error (ME) and standard deviation (Std) between measurement and fitting curve are 0.0070 and 0.0141, respectively.Furthermore, the Directive Model is fitted with measurement in various incident positions; for example, Figure 2(b) gives a scattering fitting curve when  Inci = 70 ∘ .The most suitable  and   for the scattering fitting are 0.00118 and 120, respectively.The ME and Std between fitting curve and measurement are 2.16 − 4 and 5.77 − 4, respectively.

Interpolation Algorithm in RT.
The main drawback of the RT is the high computational cost according to the complexity of 3D digital map.For dynamic channel characteristics learning, especially the HSR channel, a high time resolution (Δ =   −  −1 ) of the channel is required for the study of both large-scale and small-scale channel characteristics.In this study, the distance between TX and RX ranges from zero to one kilometer (km).As the channel at 30 GHz band is studied, the distance interval between ℎ(,   ) and ℎ(,  −1 ) should be small enough (normally less than half-wavelength of EM wave which is 5 millimeter (mm) in this study [31]).The geometry-based path interpolation is employed to overcome the impractical simulation time due to the high computational complexity [32].The basic idea of the interpolation is to obtain the information about the continuous propagation paths between two consecutive scenario snapshots.Then the linear interpolation will be performed between the two continuous paths.The detailed description of this algorithm is in [32].For this simulation in an arched tunnel, the initial time resolution is 10 milliseconds (ms) in the simulation, so the initial interval between two sampled snapshots is 1 m as the speed of HSR is assumed to be 360 km/h in this study.Then, the interpolation algorithm is preformed between each two snapshots.By utilizing geometry-based ray information, the method of interpolation will drop time resolution to 2 ms (2 mm interval in distance) in this study.Therefore, with very small time resolution, the extraction of small-scale fading parameters can be guaranteed.

Simulation Scenario and System
Setups in the Tunnel

Tunnel Scenario in Simulation.
In this study, the straight arched tunnel is employed as the HSR tunnel scenario.Figure 3 shows the overview of the tunnel scenario in the simulation.Figure 3(a) illustrated that the cross section of tunnel includes the accurate dimension and the locations of TX1 and RX1.According to realistic "Type II" tunnel described in [14], the arched tunnel in this study is with dimension 8.41 m × 6.87 m ( tunnel ×  tunnel , where  tunnel and  tunnel are defined as maximum width and height of the tunnel, resp.).The heights of TX1 and RX1 are 6.5 m and 3 m, respectively.Both TX1 and RX1 are located in the middle of the tunnel.The prior works show that totally 18 smooth surfaces which constitute the tunnel digital map can provide effective results and keep the computational complexity at reasonable level.As depicted in Figure 3(c), the train is installed with two antennas in the head and the tail, respectively; the distance between two adjacent base stations (BSs) is 1 km; the train is assumed to be 200 m long.As indicated in Figure 3(c), the tail antenna (RX1) communicates with the backward BS (TX1) while the head antenna (RX2) communicates with the forward BS (TX2).Note that, in this study, only the channel characteristics of Link 1 are investigated with the symmetric manner of Link 1 and Link 2. Furthermore, two radiation lobes of antennas in Figure 3(c) are used to illustrate the pointing directions of the directional antennas which are used in the TXs.Moreover, Figure 3(b) illustrates a snapshot in the simulation.In order to characterize the channel in this tunnel, the simulation scenario and system setups in this study follow the real requirements of the mobile communication system described in [15].The scenario and system setups of the simulation are listed in Table 1.

Antenna Setups.
As it is widely recognized that the directional antenna is indispensable for mmWave communication system, the detailed effects of directional antenna on mmWave channels still lack a careful investigation, especially for HSR channels.In other words, though the proper use of directional antenna gives a high compensation to received power, the distinctions of the radio channels characteristics between applying directional antenna and applying omnidirectional antenna are not clear.Therefore, in this study, three   antenna setups are employed under various combinations of directional antenna and omnidirectional antenna; refer to Figure 4.They are defined as follows: (1) Direc.-Direc.: TX1 and RX1 are both with the directional antennas.The antenna at TX1 is statically pointing along the tunnel while the antenna at RX1 is pointing at opposite direction of TX1 antenna.(2) Direc.-Omni.: TX1 is with directional antenna and RX1 is with omnidirectional antenna.The antenna at TX1 is statically pointing along the tunnel.(3) Omni.-Omni.: TX1 and RX1 are both with the omnidirectional antennas.
The directional antenna is designed by drawing sinusoidal chart in polar coordinate as horizontal and vertical patterns.
The following studies will be carried out in describing radio channel characteristics of three antenna setups.

Number of Frequencies in Simulation.
As the simulation configuration should follow the real mobile communication system described in [15], the simulation frequency bandwidth is set to be 125 MHz.Moreover, we want to expose the channel characteristics over a frequency range of 31.5 GHz∼33.5 GHz.So, the number of the center frequency points   is chosen sufficiently enough.Therefore, total 264 center frequency points are considered in this study [9,33,34].

Order of Multireflection in Simulation.
The computational complex of RT simulation is significantly affected by multireflection which requires massive cyclic and traversal search of RT kernel.Although higher order of reflection greatly decreases the efficiency of RT simulations, the higher order of reflection gives more accurate simulation results.Therefore, the order of reflection should be selected very carefully.Before massive RT channel simulations, we studied the power contributions of each order of reflection in the same arched tunnel in Omni.-Omni.case. Figure 5 illustrates the percentages of overall power of each order of reflection compared to total received power.In the figure, on the whole, it is apparent that the overall received power which includes power contributions from LOS to 5th order of reflection strikingly reaches 99% of total received power.Accordingly, we recorded the RT computation time of 1000 snapshots for considering up to th order of reflection (cf.Table 1).As can be seen, the computation time is of exponential growth with the increase of reflection order.Comparing results of Figure 5 and Table 1, the percentages of power contribution of reflection orders higher than 6th are less than 0.1% while they computation time are strikingly hundred times larger than that of lower reflection orders.Therefore, we limited maximum reflection order at 5th which gives accurate simulation results while keeping the computation time acceptable.
Finally, Table 2 gives an overview of scenario setups and simulation configuration in this study.

Study in System Feasibility and
Partition for Regions The value  is the received power without small-scale fading which is excluded by averaging received signal with a 40 sliding/overlapped window [35]. is the channel simulation bandwidth; it is 125 MHz in this study.  is the noise figure which is the noise factor expressed in decibel.Here,   is practically assumed to be 10 dB.Further, the number −174 is widely used as spectral noise power density for 1 Hz.As shown in Figure 6, the channel performances of different antenna setups vary considerably in the near region (which is defined in Section 4.2).In Direc.-Direc.and Direc.-Omni.cases, with the RX1 moving away from the TX1, the LOS component gradually enters the illumination of the main lobe (3 dB beam width) of the directional antenna.This process causes the SNR obvious increases in the distance ranging from 15 m to 50 m.Afterwards, when RX1 moves into far region, the LOS and lots of NLOS components enter the illumination of the main lobe of the directional antenna.
Obviously, the fluctuations of the SNRs of different antenna setups undergo a similar tendency in far region.
According to Figure 6, if the minimum SNR is 10 dB for a reliable detecting [6], that is, system minimum available threshold, this communication system in tunnel can support more than 1 km coverage range in the Direc.-Direc.and Direc.-Omni.cases.But, it is difficult to support 1 km signal coverage when deep shadow fading exists.In the Omni.-Omni.case, the system can only support coverage range less than 50 m.Although it is obvious that the directional antenna brings a better performance of signal coverage, the detailed channel characteristics are still under research.

Definition of Regions for Radio Channel Analysis.
As it is shown in Figure 6, the received SNRs are obviously different in different regions.The reasons are mainly depending on the antennas used in simulation, for example, half-power bandwidth (HPBW), pointing direction, and position [36].The following channel characteristics should be studied in different regions.Figure 7 gives the sketch of partition for regions. BS and  Ant. are the heights of BS (TX1) and RX1, respectively;  and  are the inclination angle and HPBW of the directional antenna (TX1).The red solid line indicates the pointing direction of the TX1.The two black dotted lines depict the region which will be illuminated by the antenna main lobe.The value  determines the boundary between near region and far region which can be calculated by [36] In this study,  = 8 ∘ ,  = 0 ∘ ,  BS = 6.5 m, and  Ant.= 3 m.The angle ( + /2) indicates the pointing direction plus half of the (elevation) HPBW.As a result, the length of  in this study is 50 m.

The Detailed mmWave Channel Characteristics in HSR Tunnel Scenario
Here, the radio channel characteristics will be presented in time, frequency, and polarization domains in order to help in designing a robust and sophisticated wireless communication system in HSR scenarios.

Path Loss and Shadow Fading Extraction.
The largescale fading (including path loss and shadowing fading) are obtained by averaging received power with a 40-wavelength window [35].The large-scale fading is generally expressed as a log-distance path loss model () with a path loss exponent (): where () is the function of  which indicates the distance between TX1 and RX1.( 0 ) is the intercept value at referenced distance ( 0 ).  is the shadow fading.Figure 8 gives one example of extraction process of  at a frequency center of 32.5 GHz in far region.Note that the following analyses for path loss and shadow fading are mainly in far region.This is because, in the near region, the fluctuation of the received power is largely dominated by antenna radiation pattern.
In Figure 8, a red solid line indicates the least square fitting result of the simulated data (marked in blue).The path loss exponent is the slope of the red solid line.It is around 1.1 in far region, which indicates the small attenuation of wave propagation in tunnel.This character may stem from the waveguide effects caused by tubular structure of the tunnel where the reflection attenuation will be small enough when incident angles of reflected rays are quite large in far region.The path loss exponents of different antenna setups at frequencies in the range 31.5 GHz∼33.5 GHz are calculated and shown with statistical values in Table 3.

Amplitude Distribution of Shadow Fading.
As expressed in (6), the shading fading   can be extracted from the large-scale fading ().  is conventionally modeled as a log-normal distribution [14,36] which is confirmed in this study.Figure 9 describes the probability density function (PDF) of shadowing fading in Direc.-Direc.case at 32.5 GHz with the results of  = −1.7052dB and  = 3.4443 dB, although, as the shadowing fading was extracted by (6), some deep fading inevitably leads to a no-zero value of .However, the modeled  is still valuable for studying the channel shadowing characteristics in tunnel.The statistic value of  at whole frequencies can be found in Table 3.

Autocorrelation of Shadow Fading.
As shadowing causes the channel deep fading, the communication links tend to be interrupted (refer to Figure 6).To overcome the potential communication interruption, the autocorrelation of shadow fading should be well-studied.The autocorrelation coefficient of shadow fading is one important characteristic for designing distributed antenna system, which is defined as where {⋅} denotes the expectation; () is the expression of the shadow fading at distance ; () is the expression of the standard deviation for the shadow fading at distance .Further, two widely used empirical models are employed to fit the autocorrelation coefficient [14]: the exponential model and 802.16J model.The former is accepted in WINNER II model: the latter is presented in standard IEEE 802.16J: In ( 8) and ( 9), Δ is the distance between two interested positions ( 1 and  2 ).There are mainly two definitions for decorrelation distance  cor :  cor (0.5) and  cor ( −1 ).They present the correlation coefficient equal to thresholds 0.5 and  −1 , respectively [14].Obviously, these two models have same structure.Figure 10 gives the autocorrelation coefficient of shadow fading in Direc.-Direc.case at 32.5 GHz in far region.
In conclusion, 802.16J model is fitting well when threshold is 0.5.The exponential model performs better when threshold is  −1 .

Decorrelation Distance of Shadow Fading.
Compensating the deep shadowing fading of the channel is generally used in multiantennas technology.Therefore, antennas should be separated long enough to obtain channel diversity gain.This distance is so-called decorrelation distance.
The decorrelation distances for each RX1 position along the tunnel is extracted by using (7) with two thresholds (0.5 and  −1 ) at whole frequencies in the range 31.5 GHz∼ 33.5 GHz. Figure 11 3 lists their statistic values.
For all three antenna setups, the decorrelation distances calculated by using threshold  −1 are no doubt longer than that using threshold 0.5.Moreover, Direc.-Direc.performs a longer decorrelation distance value (mean value around 2.5) than other cases (mean value around 2.3 m and 2 m, resp.).These indicate that the decorrelation distance will become longer when directional antenna is employed.

Rician 𝐾-Factor
for the Received Signal.The time-varying fading characteristic of the signal is normally modeled by Rician -factor when LOS component exists [37].the Rician -factor is defined as the ratio of the power of LOS component to the total power of NLOS components.
The Rician -factor at 32.5 GHz is shown in Figure 12(a).In the near region, the -factor experiences a rapid change.In Omni.-Omni.case, this process can be described by the narrow structure of the tunnel which causes the attenuation of reflected components.For other two cases, this rapid change mainly stems from the radiation pattern of directional antenna.In the far region, the -factors of all three antenna setups decrease slowly.Moreover, Figure 12(b) gives the CDF of three antenna setups at whole frequencies in the range 31.5 GHz∼33.5 GHz in both near region and far region.It is obvious that, in the near region, the -factor varies intensely.But, in the far region, the differences among three -factors are fairly small.The statistical values are listed in Table 3 where we find that the -factor (in dB) is positive in near region but is negative in the far region.This characteristic indicates that the dominant power contribution of the received signal is changing from the LOS component to the NLOS components.

Delay
Characteristics in the Tunnel.The root mean square (RMS) delay spread is widely known as the single parameter that can provide a quick overview of channel delay characteristics.It is defined as the normalized second-order moment of the power delay profile (PDP) which characterizes channel delay dispersion [37].In this study, the RMS delay spread is calculated as follows: where   () is the RMS delay spread;   () is the power of th ray.As all the rays are specific with certain delay, power, and angle information, ( 10) is efficient for calculating the RMS delay spread directly from rays of RT kernel results.Figure 14 depicts CDFs of the RMS delay spreads which were extracted in every snapshot at whole simulation frequencies.Apparently, when directional antennas are used in the system, the RMS delay spread will be decreased especially in near region.These results are in line with the phenomenon displayed in Figure 13 which specially compares PDPs of three antenna setups in near region.It is clear that the directional antenna can be a great spatial filter in near region that attenuates multipath components which are not illuminated by the main lobe of directional antenna.Therefore, when directional antennas are employed at both TX1 and RX1, the   minimum value of RMS delay spread can be obtained which is around 2.12 ns in near region and 0.47 in far region.

Doppler Characteristic in the Tunnel.
As discussed previously, the HSR channels in tunnel were simulated at speed of 360 km/h.Therefore, the Doppler effect on the channels is widely of interest as it gives physical interpretation of the frequency shift caused by movement [37].As can be seen in Figure 15, in near region, the train movement obviously spreads the Doppler spectrum in Omni.-Omni.case; but the spectrum shows a stable Doppler frequency shift with limited frequency spread in far region.The striking variations of the Doppler spectrum in near region are partly due to the fact that the incident angles of received rays are sparse and change rapidly, whereas, in the far region, the incident angles of received rays are very close and change slowly (referring to the tunnel narrow structure), which leads to a stable Doppler frequency.Meanwhile, since the directional antenna attenuates lots of rays in near region, Figure 16 illustrates the detailed effects of directional antenna on Doppler spectra where some distinct differences among Doppler spectra are clearly shown.To better evaluate the Doppler effects, the CDFs of mean Doppler shifts and RMS Doppler spreads of three antenna setups are studied at whole frequencies in Figure 17.These two parameters are the moments of the Doppler spectra which can be calculated similar to the moments of the PDP [37].According to Figure 17 and the statistic values listed in Table 3, the same conclusion as that from Figure 16 can be obtained that directional antenna is like a spatial filter which causes larger mean Doppler shift and lower RMS Doppler spread in HSR tunnel.
Generally, the Doppler effects are studied with other second-order fading statistics that also closely related to channel dynamic characteristics and the quality of received signals.They are level cross rate (LCR) and average fade duration (AFD) [38].As detailed in Figure 18

Polarization
Characteristics in the Tunnel.The above analyses of channel characteristics describe the vertical polarization only.However, the channel polarization characteristics in the tunnel are valuable for system design.The radiated signal from a vertically polarized antenna will experience interactions that result in horizontal polarized component before it arrives at RX and vice versa [37].This process is so-called depolarization process that has been depicted by a 2 × 2 matrix which is formulated by   ().Then, cross-polarization discriminations (XPD including XPD  and XPD  ) are invited, which indicates leakage from one polarization to another by depolarization effect [39,40]: Figure 19(a) illustrates XPD  and XPD  simulated data in Direc.-Direc.case accompanying fitting lines.The depolarization processes of the channel are found to be similar for both vertical and horizontal polarized signal at 32.5 GHz.The tubular structure of tunnel which results in the depolarization process of the channel does not obviously prefer vertical or horizontal polarized signal.Moreover, the values of XPDs in near region are much larger than values in far region.As LOS component cannot be depolarized and gradually illuminated by antenna main lobe in near region, the numerator of (11) will increase with attached antenna gain to the LOS component.This causes an apparent change as shown in Figure 19(a).
Figures 19(b) and 19(c) show the CDFs of XPD  and XPD  of three antenna setups at whole frequencies in different regions.We learn that using directional antenna can baffle depolarization process of the channel in near region.The statistic values of the results of depolarization process in the tunnel can be found in Table 3.

Relative Discussions.
Based on the radio channel characteristics and aforementioned analyses, Table 3 lists all the statistic values of the parameters.Then, the following discussions can be made: (1) Path loss: the path loss exponents for each antenna setup are all around 1.1 in far region.The tubular and narrow structure of the tunnel probably causes the waveguide propagations which compensate the received power in far region.(2) Shadow fading: as the log-normal distribution can fit the shadow fading well, the shadow fading standard deviations for three antenna setups are similar and around 3.4 dB in the far region.This number is close to the measurement results in the tunnel at 5.7 GHz [14].
(3) Decorrelation distance: the deep fading caused by shadow fading will decrease received power and decay system capacity.Also, possibly, it may result in pingpong effect when system is in handover.The decorrelation distances calculated by threshold  −1 are obviously larger than that using threshold 0.5.The mean value for Direc.-Direc.case is around 2.5 m (threshold 0.5) which is larger than other cases.In conclusion, as per directional antenna employed in system, the decorrelation distance increases 0. (5) RMS delay spread: the RMS delay spread will be deceased when directional antenna is employed.In the near region, the directional antenna performs as a spatial filter which reduces the RMS delay spread by 1-1.5 ns.In the far region, RMS delay spreads of different antenna setups are less than 1 ns, and in the Direc.-Direc.case, the number will be less than 0.5 ns.
(6) RMS Doppler spread: there are less than 100 Hz variations in RMS Doppler spreads when one or two directional antennas are employed in the system.All RMS Doppler spreads in far region are around 0.41 kHz.Furthermore, the RMS Doppler spread may be less influenced than RMS delay spread when directional antennas are employed.case, the values of XPD  are 2 dB larger than XPD  in both the two regions.Since the vertical polarized signal shows an advantage in resisting depolarization compared to horizontal polarized signal in the tunnel, it is better to deploy vertical polarized antennas in the arched tunnel.

Conclusion
In this paper, 30 GHz band HSR radio channels in an arched tunnel with different antenna setups are investigated by a RT tool.To enhance RT for more accurate simulations, a special material measurement campaign is performed for the cement which is the main construction material of the tunnel.Then, an advanced time-interpolation method is described for extracting small-scale channel characteristics.These two steps work as a promising and solid base for our extensive channel simulation at 31.5 GHz∼33.5 GHz.The channel characteristics are found to be widely different in different regions.The results show that directional antennas deployed at both TX and RX can significantly improve the coverage range of mobile communication system.But the directional antennas will bring about larger decorrelation distance in radio channel.Therefore, the interval between antennas at RX (if multiantennas are applied) should be separated by more than 2.5 m to overcome deep shadow fading.Meanwhile, the directional antennas play the role of spatial filter which obviously attenuates the rays out of the antenna main lobe, especially when RX is close to TX.These effects greatly decrease RMS delay spread in the near region.Actually, comparing to omnidirectional antenna, the directional antenna mostly affects the channel characteristics in the near region, but it rarely affects the channel characteristics in the far region.Considering the length of the near region is around 50 m in the study, the HSR will get through the near region in a flash (0.5 second).In this moment when the channel changes with extreme rapidity, it is difficult to keep a reliable communication link.So, advanced handover strategies should be considered before or after HSR getting through the near region.Finally, the vertical polarized antenna is suggested in this arched tunnel rather than horizontal polarized antenna.

Figure 1 :
Figure 1: Study on EM and scattering parameters of cement wall.(a) Geometry of material study where the tunnel image is from [30].(b) Measurement campaign.

Figure 2 :
Figure 2: The comparison between fitting curve and measurement of cement brick at 32.5 GHz.(a) Reflection.(b) Scattering.

Figure 3 :
Figure 3: (a) Cross section of the tunnel in the simulation.(b) One snapshot of the simulation in the tunnel.(c) The sketch map for the scenario.

Figure 6 :
Figure 6: SNRs of different antenna setups.The red vertical solid line is the partition for two regions; the black horizontal line depicts the minimum threshold for the available SNR.

Figure 8 :
Figure 8: Path loss fitting at 32.5 GHz in Direc.-Direc.case in far region.

Figure 9 :
Figure 9: Fitting for amplitude distribution of shadow fading.

Figure 11 :
Figure 11: (a) Decorrelation distances of three antenna setups at 32.5 GHz in far region.(b) For three antenna setups, the figure shows the CDF of the decorrelation distances at whole frequencies with thresholds 0.5 and  −1 in far region.

Figure 12 :
Figure 12: (a) Rician -factor for different antenna setups at 32.5 GHz.(b) CDF of Rician -factor in different antenna setups and regions at whole frequencies.

Figure 13 :
Figure 13: Partial enlarged views to describe the changing of the PDPs with different antenna setups in near region at 32.5 GHz.

Figure 14 :Figure 15 :
Figure 14: CDF of RMS delay in different antenna setups and regions at whole frequencies.
(a), curves that stand for LCRs in far region exhibit lower LCRs compared to that in near region; and the result of Direc.-Direc.case in far region shows the minimum LCR.These phenomena are in line with observation results in Doppler spectra (cf.Figures 16 and 17).Additionally, since the LCRs of far region are far lower than LCRs of near region, the Doppler frequencies in far region are accordingly more stable with minimum Doppler spread and mean shift.However, although the AFDs exhibited in Figure 18(b) have no clear distinctions between different regions, the AFD of Direc.-Direc.case still shows the minimum value.And AFDs of three antenna setups in far region increase sharply when thresholds for AFD calculation are around 0 dB, which indicate weaker fastfading characteristics of the channels in far region as well as weaker Doppler spread.

Figure 16 :Figure 17 :
Figure 16: Partial enlarged views to describe the changing of the Doppler spectra with different antenna setups in near region at 32.5 GHz.

2 - 0 .4 m. ( 4 )
Rician -factor: the -factors in two regions show great differences, especially for Direc.-Direc.case where the difference between the mean values of two regions is around 25 dB.Even if the -factors have a great variation in different regions, the -factors of different antenna setups have little differences in far region (around −6 dB).

( 7 )
XPDs: for Direc.-Direc.and Direc.-Omni.cases, it is interesting to find that the mean values of XPD  and XPD  are very close in near region and have a variation of 1.7 dB in far region.For Omni.-Omni.

Table 1 :
Computation time (min) of total 1000 snapshots for orders of reflection.

Table 2 :
Scenario setups and simulation.