#### Abstract

This paper presents a three-dimensional (3D) geometry-based stochastic model (GBSM) for capturing the non-stationarity of radio channel at 1.8 GHz in a rectangular tunnel. A time-variant (TV) complex channel gain is derived for obtaining the statistical properties in time, frequency, and spatial domains such as the time-variant autocorrelation function (TV-ACF), the time-variant Doppler power spectral density (TV-DPSD), and the time-variant spatial cross-correlation function (TV-CCF), respectively. Then the TV channel statistical properties at different time instants can be extracted and the non-stationary channel characteristics caused by the TV scattering environment are thoroughly discussed. Furthermore, three cases including “approach”, “arrival”, and “away” are set to allow a comprehensive study on how the DPSD behaves with the relative position between transmitter and receiver. The reliability of proposed 3D GBSM is highlighted by a good agreement with the measured result in terms of the correlation function.

#### 1. Introduction

For realizing a more stable and reliable train-ground communication, the multiple-input multiple-output (MIMO) technology working at 1.8 GHz has been gradually utilized in metro communication system. The MIMO channel model is thus a prerequisite for predicting the coverage of wireless link and optimizing the system performance. The geometry-based stochastic model (GBSM) [1–3] which combines the geometry of the environment and the statistical information is regarded as a flexible way to describe the MIMO channel characteristics. Considering the propagation in the elevation plane, lots of attention has been paid to three-dimensional (3D) channel models [4–6]. In [4], a concentric-cylinder model was derived for wideband MIMO mobile-to-mobile channel, and the independent azimuth angle and elevation angle was analyzed separately. To study the joint effect of the azimuth and elevation angle, 3D two-sphere models have been proposed in [5, 6]. Several GBSM based standard channel models also have been proposed, e.g., 3GPP-3D model, IMT-Advanced model, and WINNER+ model. These standard models were built with the angle parameter obeying Laplacian, Gaussian, and Wrapped Gaussian distributions, respectively. Furthermore, these models were evolved to support to non-stationary scenario by replacing the model parameters with the time-variant (TV) ones [7, 8]. In spite of this, these standard models are still focused on propagation scenarios like Indoor, Urban, Rural, and Suburban; mobile propagation scenarios are not involved in these models. Another modeling approach in [9, 10] employed random mobility trajectories of the mobile station (MS) to inject the non-stationarity into the characteristics of the channel, but without the consideration of the regular movement in confined space, especially in tunnel scenario. It is worth mentioning that the radio signal propagations inside tunnel will suffer more reflections and scattering due to the confined space, rough interior tunnel walls, and waveguide effect [11]. To acquire better performance of metro communication system, it is essential to build accurate tunnel channel models to comprehensively understand the statistical properties of tunnel channels.

Few studies were devoted to the non-stationary channel models for the tunnel environment. The regular-shaped GBSM in cutting and irregular-shaped GBSM in suburb and urban scenarios were proposed in [12]. In [13], a geometry-based single-bounce model of MIMO channel in the tunnel has been proposed in order to investigate the factors which influence the channel correlation, but it is estimated without considering the elevation angles. A geometry-based MIMO car-to-car channel model is applied for a semi-circular tunnel [14]. However, this model only considered the local scatterers between two cars which both move at the same speed. The local angle of arrival (AOA) and angle of departure (AOD) were thus time-invariant, and the non-stationarity was unable to be characterized. The two-dimensional (2D) high-speed train model in [15] considered the non-stationary characteristic caused by the variation of the velocity and the angle of motion (AOM), but the temporal variations of AOA which contribute to the non-stationary characteristic were neglected.

However, the modeling of AOA and AOD plays an import role in the non-stationary mobility model. Empirical stochastic distributions, e.g., the uniform distribution [16, 17], Von Mises distribution [18, 19], and Von Mises-Fisher distribution [20], have always been assumed to characterize the angular parameter to achieve the closed-form representation of channel statistical properties. But the AOAs in [16–20] are not modeled as time-variant functions. It is not suitable for the real-time mobile behavior modeling; the temporal variations of the angle parameters need to be included in these models.

Inspired by the issues mentioned before, we present a 3D geometrical model to describe the non-stationary channel caused by the movements of the MS in the tunnel environment. By introducing the TV AOA and propagation length, we derive the TV channel gain based on the single-bounce scattering. Furthermore, in the non-stationary process, the phase of the complex channel gain is no longer rotated with the TV Doppler frequency multiplied by time, but the integral of the TV Doppler frequency. The former multiplicative relationship can only be obtained by developing the TV phase in a second-order Taylor series around* t* = 0 [21]. By employing the analytical method of the non-stationary processes, the expressions and analytical results are given for the time-variant autocorrelation function (TV-ACF), the time-variant Doppler power spectral density (TV-DPSD), and the time-variant spatial cross-correlation function (TV-CCF). Moreover, the correctness of the derivations is verified by the simulation results. Finally, a measurement is conducted in a subway-like tunnel and the measured results are used to validate the effectiveness of the proposed model.

The paper is organized as follows. Section 2 proposes a new 3D non-stationary GBSM, and the statistical properties of the model are derived. The theoretical results are depicted and compared with the corresponding simulation in Section 3. In Section 4, the measurement results are discussed with the proposed model. Finally, Section 5 concludes the paper and outlines future research lines.

#### 2. The 3D Non-Stationary GBSM

A typical propagation scenario for the metro environment is shown in Figure 1. The downlink is from the base station (BS) to the moving MS. The symbols of and denote the BS transmitter (Tx) and the MS receiver (Rx), respectively. The proposed model generates the empty straight tunnel scattering environment with a length of . The tunnel cross section is rectangular with a width of and a height of . Figure 2 shows the distribution of scatterers on the inner surfaces of the tunnel wall. In order to study the general cases of the scattering inside the tunnel with the moving mobile station, it is assumed that the material of the tunnel wall is uniform and the scatterer distribution follows uniform distribution.

The 3D MIMO non-stationary GBSM is described as follows. The measurement inside the tunnel [22] demonstrates that most of the propagation wave energy was from the line-of-sight (LOS) and single-bounced (SB) components. Therefore, the LOS and SB components are under consideration in our proposed model. The signals received at the moving Rx are reflected from the infinite scatterers which are randomly distributed along the walls, ceiling, and ground of tunnel. For the simplified description, (, , ) are used to describe the positions of the scatterers on the inner surfaces of the tunnel wall for , , and . Since the upper limit* M*,* N*,* O* with each tending to infinity, three discrete random variables , , and can be denoted as the continuous variables* x*,* y,* and to acquire the theoretical model.

##### 2.1. Channel Gain

The channel gain between the* p*-th Tx element and the* q*-th Rx element is demonstrated as the superposition of the LOS and the SB components, and it can be expressed as

For the SB component, the channel gain can be formulated as the summation of several multipath components from scatterers as given belowwhereThe key parameters in equations are defined in Table 1. The symbol in (2) represents the Rice factor, which is defined as . denotes the mean power between the* p*-th Tx element and the* q*-th Rx element. In (4), and are used to represent the positions of the and . Using the trigonometric relationship, the azimuth angle of arrival and elevation angle of arrival can be expressed asThe azimuth angle of departure and the elevation angle of departure can be derived in the similar way. In (6), indicates the maximum Doppler shift and denotes the TV Doppler shift of the SB component caused by the movement of the Rx, where is the wavelength.

The LOS component of the total channel gain in (1) can be written as whereIt is worth mentioning that the in (2) and the in (9) are the initial phases of the SB and the LOS component; each is assumed with uniform distribution over the interval . The symbol in (11) denotes the TV Doppler shift of the LOS component caused by the movement of the Rx. The TV angle parameters , , and also can be computed in the similar way as (7) and (8).

##### 2.2. Statistical Properties of the Model

The typical TV channel statistical properties can be derived from the channel gain, such as the TV-ACF, TV-DPSD, and TV-CCF.

###### 2.2.1. TV-ACF

The TV-ACF of the subchannel from the* p*-th Tx element to the* q*-th Rx element is defined bywhere ∗ denotes the complex conjugate operator and is the statistical expectation operator.

Substituting (2) and (9) into (12) and assuming that the propagation distance is invariant for short time difference and the phase rotated with the integral of the TV Doppler frequency has a prominent impact on the TV-ACF, the TV-ACF of the LOS and the SB component can be obtained as follows:where is the joint probability density function (PDF) of the continuous random values* x*,* y*, and* z*.

###### 2.2.2. TV-DPSD

For non-stationary channels, the TV-DPSD can be obtained from the TV-ACF by applying the Fourier transform with regard to time difference , as given below

It can be executed to reflect the variations in the power spectrum along the Doppler frequency throughout the movement of the Rx.

###### 2.2.3. TV-CCF

The normalized TV-CCF between two arbitrary channels in the spatial domain isFor the case of the LOS component, it can be defined asFor the case of the SB component, it can be defined asMoreover, the CCF at and can be denoted as and , respectively. As an example, by substituting in (18), can be obtained as follows:

The uniform distribution is assumed in scattering model; therefore the joint PDF can be computed as [23]

#### 3. Simulation Results and Validation

Here, we investigate the statistical properties of the proposed model in detail. A 150 m long rectangular tunnel with the height of 5 m and the width of 5.2 m is considered.

##### 3.1. Simulation Model

By assuming the infinite number of the scatterers (), the analytical model described in (2) is hard to realize. However, a corresponding simulation model with a finite number of scatterers can be developed. In terms of the model parameters, the discrete random values are generated by discrete uniform random parameters, while the other parameters are the same as the theoretical model, summarized in Table 2. To achieve a better approximation to the statistical properties of the theoretical model, the simulation results are obtained by averaging the results of multiple simulation runs (up to 10^{4}), and the number of scatterers per simulation run is selected as . Furthermore, the statistical properties of the simulation model are shown as (21)-(23). For TV-ACF, the simulation model is expressed asFor TV-DPSD, the simulation model is given asFor TV-CCF, i.e., the CCF at , the simulation model is written as follows:

##### 3.2. TV-ACF

By using (12)-(14), the numerical results obtained for the TV-ACF in NLOS scenario are shown in Figure 3. For clearly comparison purpose, the absolute values of the theoretical TV-ACF at different time instants* t* = 0 s, 1.5 s, 3 s (corresponding to = 0 m, 33.3 m, 66.6 m) are extracted as shown in Figure 4. Also, the corresponding simulation results can be obtained by (21). It can be seen that a good agreement between the theoretical and simulated results is achieved, which verifies the derivation of the proposed model. Furthermore, different variation trends of the ACF are resulted from the dynamic instantaneous phase which varies with the TV Doppler frequency. As the Rx gets closer to the Tx (fixed location, = 75 m), the channel autocorrelation behaves more sensitive to the time difference (see, for example, the curve at* t* = 3 s). From Figure 5, it is obvious that the TV-ACF of the channel with the LOS component is higher than the case without the LOS component. In addition, the non-stationarity of channel characteristic is not notable in LOS dominated propagation condition.

##### 3.3. TV-DPSD

The theoretical TV-DPSD calculated by (15) is given in Figure 6. With the variation of the angle between the incoming wave and the direction of motion, the Doppler shift has a smooth drift from positive maximum to negative minimum. Figure 7 depicts three typical Rx locations with respect to Tx, the approach location (*t* = 0.45 s, = 10 m, Case I), the arrival location (*t* = 3.375 s, = 75 m, Case II), and the away location (*t* = 6.30 s, = 140 m, Case III). They are extracted from Figure 6 for further illustration as shown in Figure 8. A symmetric DPSD can be obtained in Case II while asymmetric DPSDs exist in Case I and Case III. In the non-isotropic environment like Case I and Case III, the AOA has a relatively concentrated distribution only in one direction. As a result, the energy focuses on the positive maximum Doppler shift and negative minimum Doppler shift, respectively. However, the reason why the energy is approximately uniform distributed in Case II is that the propagation scenario can be assumed to be a rich scattering environment in which the directions of incoming waves are uniform. In addition, the theoretical results fit well with the corresponding simulation ones obtained by (22).

##### 3.4. TV-CCF

The absolute value of the TV-CCF can be computed using (16)-(19), as depicted in Figure 9. It can be easily observed that the TV-CCFs decrease as the antenna element spacing of the Rx increases. The TV-CCFs at different time instants* t* = 0 s, 1.5 s, 3 s (corresponding to = 0 m, 33.3 m, 66.6 m) are extracted in Figure 10, respectively. It is demonstrated that when the Rx gets closer to the Tx (fixed location, = 75 m), the TV-CCF decreases more sharply with the increasing of the antenna spacing. This trend is caused by the variation of multipath dispersion in the spatial domain. When the Rx is far from the Tx, the multipath is relatively concentrated in one direction, while it is getting more dispersed as the Rx closing to Tx. A lower CCF is friendly to the MIMO system working in the way of spatial diversity in which the sub-channels are expected to keep independent with each other. Moreover, the corresponding simulation results are computed utilizing (23) and the good correspondence between the theoretical and simulated results validates the correctness of the above derivations.

#### 4. Experimental Verification with Tunnel-Like Environment Measurement

A measurement campaign was carried out in a tunnel-like scenario at Zhongtian Technology Company, Nantong, China. As shown in Figure 11, the overall length of the experimental tunnel is 100 m, with 50 m long rectangular part and 50 m long arched part. A concrete section exists in the junction which combines these two parts as shown in Figure 11(a). The 50 m long rectangular part has been considered for the measurement. The size of the rectangular cross section is 4.4 m (width)×3 m (height). The measurement was performed in three experiment regions denoted as A1, A2, and A3. In each area, a matrix of 1 × 8 receiving points was configured to simulate a virtual receiving array, as shown in Figure 11(b). The transmitting antenna was set at the middle of the tunnel, which was 0.5 m far away from the tunnel wall and 1.6 m high as illustrated in Figure 11(c). In the processing of the measured data, the channel impulse response (CIR) is directly generated by the cross-correlation of the measured data collected by the Rx with the replica of the transmitted PN sequence. Then the multipath components are detected with significant power contributions above the 5 dB of the estimated multipath noise floor. And finally, the narrowband CIR obtained by the complex sum over the delay domain is substituted into the expression (10) in [24] to extract the spatial correlation coefficient. As for the theoretical model, the spatial CCF can be obtained by (23) and the model parameters are set according to the measurement campaigns. Figure 12 shows that the theoretical curves can properly capture the measurement results for spatial CCF in A1, A2, and A3; especially a good performance is found in A2. This result demonstrates that the proposed model can generate a similar scattering environment as a realistic tunnel-like environment.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

#### 5. Conclusion

In this paper, a 3D GBSM for non-stationary channels in tunnel environment has been proposed. By introducing the time-varying multipath length and AOA, the proposed model can be applied to simulate the non-stationary channels in real time. The different variation trends of TV-ACF at different time instants are related to the dynamic changes of the instantaneous phase in temporal. And it is obvious that the absolute values of TV-ACFs are higher in the propagation environment with the LOS components. Additionally, the TV-DPSDs for three typical cases have shown that the fluctuating trend of energy distribution in Doppler frequency arises from the variations of the scattering environment. It is also demonstrated that the TV-CCFs decrease with the increment in time, and when the multipath in the spatial domain is relatively dispersed, the decreasing rate is higher. Finally, a good fitting of the spatial CCFs of the theoretical model with the measured results verifies the feasibility of the proposed model.

In future studies, the moving scatterers or the variation of the scatterers' distribution should be incorporated in the non-stationary model. For the dedicated cases in other tunnel scenarios, we will be devoted to establishing the applicable models to construct the spatial distribution of scatterers. Validation of the proposed channel model with other real-time measurement scenarios and with respect to other channel statistical properties is also an important topic for future research. Furthermore, the proposed channel model can also be studied under other signal frequency, tunnel curvature, and moving trajectory.

#### Data Availability

The measurement data used to support the findings of this study have not been made available because of the protection of intellectual property rights. Moreover, the data is shared with our partner company and will be used in our future work; also the measurement work costs a lot of man power and equipment costs.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grants 61571282 and 61871261 and the Open Project Program of the State Key Laboratory of Rail Traffic Control and Safety under Grant RCS2017K012.